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expensive,  inefficient  trial-and-error  search.  In  order  to  expedite  such  efforts  we  have  developed  and  applied  computational  quantitative 
structure-activity  relationship  models  to  the  prediction  of  AChE-binding  efficacy  and  toxicity  for  candidate  prophylactics,  identifying  six 
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Reliable  Prescreening  of  Candidate  Nerve  Agent  Prophylaxes  via  3D  QSAR 

Final  Report 

Gerald  H.  Lushington,  Nora  M.  Wallace  and  Jian-Xin  Guo 
The  University  of  Kansas 

Forward:  This  work  was  aimed  at  developing  and  implementing  a  reliable,  efficient 
computational  method  for  predicting  the  propensity  of  specific  chemicals  to  serve  as 
prophylaxes  against  nerve  agent  poisoning.  Based  on  three  dimensional  (3D) 
quantitative  structure-activity  relationship  (QSAR)  criteria,  the  method  was  designed  to 
identify  those  Acetylcholinesterase  (AChE)  inhibitors  that  either  bind  too  weakly  (i.e.,  are 
ineffectual)  or  too  strongly  (i.e.,  are  themselves  toxic)  from  the  myriad  of  known  AChE 
inhibitors.  The  resulting  tool  should  thus  expedite  DOD  prophylaxis  research  by 
screening  out  compounds  with  inappropropriate  inhibition  tendencies. 

This  project  aims  to  deliver  software  capable  of  deriving  quantitatively  reasonable 
binding  affinity  predictions  from  docked  ligand-receptor  complexes,  thus  allowing 
estimation  of  potential  appliation  to  nerve  agent  prophylaxis.  Conventional  docking 
scoring  functions  whose  entropic  term  has  been  trained  across  a  broad  range  of  enzyme 
receptors  seem  to  score  AChE  interactions  very  poorly.  To  properly  account  for  AChE 
entropy,  we  will  use  QSAR  techniques  to  train  an  AChE-specialized  scoring  model  as  a 
weighted  linear  sum  of  ligand  atom  -  receptor  residue  enthalpic  interactions.  The 
predictions  generated  by  this  software  need  not  replace  more  rigorous  experimental  or 
computational  measurements,  but  should  permit  prescreening  that  will  help  researchers 
to  focus  their  efforts  on  manageable  subsets  of  the  extensive  collection  of  known  or 
suspected  AChE  inhibitors.  This  work  is  also  expected  to  provide  an  efficient  means  for 
intuiting  new  inhibitor  scaffolds  by  elucidating  the  steric  and  electrostatic  importance  of 
specific  receptor  residues  whose  interactions  must  be  accomodated  by  bound  ligands. 

Tables  and  Figures: 

-  Figure  1 .  Schematic  showing  the  information  flow  into  and  within  our  AChE  ligand 
binding  affinity  prediction  (scoring)  software. 

-  Table  1 .  Predicted  AChE  binding  affinity  for  71  molecules  from  the  PubChem  Small 
Molecules  Collection  with  physicochemical  attributes  similar  to  known  AChE 
inhibitors. 

-  Figure  2.  Top  scoring  AChE  inhibitors  from  derived  from  PubChem  analogs  of 
known  anticholinesterases 

-  Figure  3.  Bound  conformations  of  molecules  1 144141  (A),  and  101917  (B)  within  the 
AChE  active  site. 

Problem  Statement:  identification  of  strong  but  non-toxic  AChE  inhibitors  is  critical  to 
continued  efforts  to  devise  a  safe  and  reliable  prophylaxis  scheme  for  countering  the 
toxic  potential  of  nerve  agent  chemical  weaponry.  In  vitro  and  in  vivo  pursuit  of  such 
research  is  inherently  dangerous,  time  consuming  and  expensive,  thus  the  development 
of  reliable  new  computational  methods  to  prescreen  potential  inhibitor  candidates  is  of 
great  potential  value.  Given  multiple  crystal  structures  of  AChE  available  in  the  literature 
and  within  the  protein  databank,  this  research  is  highly  amenable  to  structure-based 
design  techniques,  however  the  inhibition  scoring  models  provided  with  conventional 
molecular  docking  software  do  not  correlate  well  with  experimental  affinity  trends  among 
known  AChE  inhibitors,  Likely  due  to  an  inaccurate  or  inadequate  account  of  entropic 
effects  within  the  AChE  receptor.  We  have  thus  a  scoring  model  in  this  work  to  reliably 


account  for  such  effects  and  thereby  reproduce  known  behavior  and  predict  efficacy  of 
hitherto  uncharacterized  compounds. 


Important  Results:  The  basic  methodology  underlying  this  work  was  developed, 
validated  and  published  [1]  prior  to  commencement  of  this  project.  The  corresponding 
publication  is  provided  as  Appendix  A.  The  primary  deliverable  in  this  project,  the 
development  of  software  implementing  the  resulting  AChE-inhibition  specific  score 
function  for  use  by  chemical  defense  researchers  (described  below)  was  completed  by 
9/30/2005.  Use  of  the  resulting  software  entails  a  process  flow  as  depicted  in  Figure  1 , 
whereby  a  researcher  applies  our  specially  tailored  AChE-specific  score  function  to 
evaluate  the  inhibitive  potential  for  ligands  that  have  been  computationally  or  manually 
docked  to  the  AChE  receptor.  The  actual  computational  docking  simulation  is  a 
prerequisite  to  our  analysis  and  is  not  explicitly  performed  via  the  software  developed 
herein,  although  our  method  is  compatible  with  the  docking  results  attained  by  most 
available  docking  programs,  including  commercial  codes  such  as  FlexX  and  MOE-Dock, 
plus  academic  freeware  codes  such  as  AutoDock.  An  additional  co-requisite  step,  the 
assignment  of  Merck  Molecular  Force  Field  (MMFF94)  atom  types  and  partial  charges 
[2]  must  also  be  performed  independently  of  our  software,  as  can  be  readily 
accomplished  by  a  variety  of  3rd  party  programs  such  as  QUACPAC,  SYBYL,  MOE,  etc. 
Given  the  docked  structure  and  the  proper  MMFF94  atom  types,  one  may  then  apply  our 
new  programs  compute  a  pKj  score  for  the  ligand  via  the  expression  described  in  detail 


Fig.  1.  Schematic  showing  the  information  flow 
into  and  within  our  AChE  ligand  binding  affinity 
prediction  (scoring)  software. 


Compatible  Software 


Ligand  model  preparation: _ 

Chem3D  (commercial) 
http://www.cambridgesoft.com 
ChemSketch  (free  for  non-profit) 
http://www.acdlabs.com 
Insight-ll  (commercial) 
http://www.accelrys.com 
MOE  (commercial) 
http://www.chemcomp.com 
MOLDA  (free  for  non-profit) 
http://www.molda.org 
SYBYL  (commercial) 
http://www.tripos.com 

Docking: _ 

AutoDock  (free  for  non-profit) 
http://www.scripps.edu/mb/olson 
FlexX  (commercial) 
http://www.tripos.com 
MOE-Dock  (commercial) 
http://www.chemcomp.com 

Atom  typing  /  charge  calculation: 
QUACPAC  (free  for  non-profit) 
http://www.eyesopen.com 
SYBYL  (commercial) 
http://www.tripos.com 


in  Eq.  1  of  Appendix  A.  Specifically,  ligand  /  receptor  residue  electrostatic  and  van  der 
Waals  intermolecular  interaction  enthalpy  terms  are  computed  according  to  the  MMFF94 
force  field  terms  as  devised  by  Thomas  Halgren  [1],  These  enthalpies  are  translated  to 
final  predicted  pKj  affinity  predictions  by  incorporating  entropic  perturbative  corrections 
(for  each  receptor  residue  within  a  12  A  sphere  of  the  central  cavity)  as  evaluated  by 
partial  least  squares  (PLS)  fitting  to  affinity  trends  within  a  set  of  53  known  non-covalent 
AChE  inhibitors.  In  carrying  out  this  PLS  modeling,  we  achieved  an  excellent  correlation 
(R2  =  0.89)  relative  to  experimental  affinity  data,  as  well  as  strong  predictivity  (leave-one- 
out  cross  validated  correlation  of  Q2  =  0.71  within  the  training  set  and  a  predictive  R2  = 
0.69  value  for  a  set  of  16  compounds  left  out  of  the  original  training  set). 

Work  during  the  final  three  months  of  CY2005  was  devoted  to  the  application  of 
our  method  to  elucidation  of  new  potential  inhibitors.  To  accomplish  this,  then  entire 
PubChem  small  molecule  virtual  library  (millions  of  compounds)  was  scanned  via  the 
DiverseSolutions™  [3]  program  to  elucidate  recently  synthesized  compounds  with 
physicochemical  properties  similar  (but  not  identical)  to  those  known  AChE  inhibitors  in 
the  original  69  compound  training  and  test  sets  reported  in  the  previous  paragraph. 
From  this  analysis,  a  total  of  71  potential  new  inhibitor  candidates  were  elucidated 
according  to  minimal  physicochemical  distances  as  scored  by  BCUT  diversity  metrics 
[4].  This  set  of  compounds  were  then  docked  into  an  AChE  receptor  model  via  the 
FlexX  program  [5],  and  computational  pKi  predictions  were  then  performed.  The 
resulting  affinity  scores  are  reported  in  Table  1,  and  the  five  top  scoring  (sub-nanomolar) 
inhibitors  are  depicted  in  Fig.  2. 


Table  1.  Predicted  AChE  binding  affinity  for  71  molecules  from  the  PubChem  Small 
Molecules  Collection  with  physicochemical  attributes  similar  to  known  AChE  inhibitors, 
potential  sub-nanomolar  inhibitors  (pKi  >  9.0)  are  shown  in  bold  font. 

Cpd.  ID  pKi  Cpd.  ID  pKj  Cpd.  ID  pK{  Cpd.  ID  pK{ 


3389 

138357 

184705 

191042 

210342 

249033 

266808 

270535 

281236 

376530 

376534 

409286 

427257 

434434 

493838 

495057 

498062 

625525 


6.4707 

7.4426 

8.9626 

7.6288 

8.7880 

10.8281 

4.3691 

5.3963 

7.0434 

8.6967 

8.8889 

8.1989 

3.5223 

6.2186 

6.9859 

6.9699 

7.4259 

8.7682 


680250 

720003 

748308 

748321 

775240 

783427 

783500 

785217 

798534 

802115 

802520 

833602 

838412 

848206 

858800 

881012 

899731 

919226 


8.6312 
5.6345 
10.1534 
5.7561 
3.3730 
7.0665 
7.7037 
7.4220 
7.4776 
7.8668 
1 .2580 
6.2348 
7.7260 
7.1518 
8.9556 
6.8794 
8.7880 
7.5994 


933051 

951510 

956831 

956837 

986516 

1001917 

1133586 

1144141 

1150763 

1152275 

1152327 

1154309 

1163388 

1167962 

1179804 

1236539 

1316823 

1333763 


5.7801 

6.9074 

5.5098 

6.4678 

4.5478 

10.6315 

7.9118 

14.2625 

7.5349 

7.2546 

7.2546 

6.0770 

6.1331 

6.7236 

6.5368 

6.8997 

7.0958 

6.2979 


1339872 

1341521 

1379265 

1410811 

1427301 

1430079 

1464709 

1467950 

1563954 

1586860 

1629088 

1637355 

1637358 

1637361 

1639895 

1658324 

1665392 


7.3119 

7.6272 

7.2713 

5.4950 

8.1817 

5.7702 

6.9736 

-8.0860 

6.7446 

5.9723 

8.4024 

6.6228 

7.3683 

7.2771 

9.0748 

7.2574 

8.4541 


Among  the  top  scoring 
prophylactic  candidates  that  we  identified, 
the  one  predicted  to  be  the  most  powerful 
inhibitor,  molecule  #1144141,  is  a  tacrine¬ 
line  peripheral-active-site  (PAS)  binder 
whose  docked  conformation  is  reported  in 
Fig.  3a.  Molecules  249033  and  748308, 
and  1639895  have  similar  PAS-binding 
conformations.  Molecule  #101917,  on  the 
other  hand,  binds  one  tail  in  the  PAS 
region  and  another  down  near  the  base  of 
the  AChE  gorge  in  a  manner,  shown  in 
Fig.  3b,  that  is  analogous  to  the  well 
known  Alzheimer's  disease  drug  aricept 
(donepezil). 


Fig. 2.  Top  scoring  AChE  inhibitors  from 
derived  from  PubChem  analogs  of  known 
anticholinesterases. 


\ 


Fig.  3.  Bound  conformations  of  molecules  1144141  (A),  and  101917 
(B)  within  the  AChE  active  site.  The  receptor  surface  is  colored  as 
follows:  red  =  H-bond  acceptors  (O's  and  lone  pair  N's),  blue  = 
donatable  H's  (HO  and  HN),  white  =  other  polar  atoms,  yellow  = 
nonpolar  atoms.  Ligands  are  colored  according  to  standard  CPK 
coloring  (H  =  cyan,  C  =  white,  N  =  blue,  O  =  red,  S  =  yellow). 

Although  molecule  #1144141  is  clearly  predicted  to  be  the  stronger  inhibitor,  it 
has  two  potential  drawbacks:  a)  while  it's  predicted  affinity  may  be  a  significant  over¬ 
estimation  (since  there  were  no  inhibitors  in  the  training  set  within  orders  of  magnitude  of 
the  same  affinity)  it  may  still  be  too  powerful  an  inhibitor  and  may  thus  effect  an 
intolerable  level  of  AChE-inhibition-related  toxicity,  and  b)  deriving  most  of  its  inhibitive 
potential  from  PAS  hydrophobic  interactions,  like  tacrine,  may  significantly  diminish  the 
oral  availability  and  distributive  potential  of  the  molecule.  Consequently,  molecule 
#101917  and  its  analogs,  essentially  exhibiting  binding  conformers  similar  to  aricept  but 
with  better  hydrogen  bonding  interactions  with  the  AChE  receptor,  may  prove  to  be  more 
fruitful  in  yielding  practical  prophylactic  prospects. 


Finally,  an  important  complementary  method  has  been  developed  to  recognize 
and  predict  prospective  toxicity  among  covalent-binding  AChE  inhibitors  of  potential 
application  to  nerve  agent  prophylaxis  and  therapeutics.  This  work  has  recently  been 
published  in  the  journal  Chemical  Research  on  Toxicology,  and  the  corresponding  paper 
is  included  herein  as  Appendix  B.  While  the  initial  work  in  this  field  focused  on  covalent 
AChE  inhibitors  due  to  a  greater  availability  of  in  vivo  toxicity  data,  modifications  to  the 
methodology  are  possible  to  extend  its  scope  to  non-covalent  systems  such  as  those 
reported  above. 
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A  Docking  Score  Function  for  Estimating  Ligand- Protein  Interactions: 
Application  to  Acetylcholinesterase  Inhibition 
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Acetylcholinesterase  (AChE)  inhibition  is  an  important  research  topic  because  of  its  wide  range 
of  associated  health  implications.  A  receptor-specific  scoring  function  was  developed  herein 
for  predicting  binding  affinities  for  human  AChE  (huAChE)  inhibitors.  This  method  entails  a 
statistically  trained  weighted  sum  of  electrostatic  and  van  der  Waals  (VDW)  interactions 
between  ligands  and  the  receptor  residues.  Within  the  53  ligand  training  set,  a  strong  correlation 
was  found  (R2  =  0.89)  between  computed  and  experimental  inhibition  constants.  Leave-one- 
out  cross-validation  indicated  high  predictive  power  ( Q 2  =  0.72),  and  analysis  of  a  separate 
16-compound  test  set  also  produced  very  good  correlation  with  experiment  (R2  =  0.69).  Scoring 
function  analysis  has  permitted  identification  and  characterization  of  important  ligand-receptor 
interactions,  producing  a  list  of  those  residues  making  the  most  important  electrostatic  and 
VDW  contributions  within  the  main  active  site,  gorge  area,  acyl  binding  pocket,  and  periferal 
site.  These  analyses  are  consistent  with  X-ray  crystallographic  and  site-directed  mutagenesis 
studies. 


Introduction 

Acetylcholinesterase  (AChE)  is  an  enzyme  that  hy¬ 
drolyzes  the  neurotransmitter  acetylcholine  (ACh)  at 
cholinergic  synapses,  accomplishing  its  role  at  a  rate 
faster  than  those  of  most  other  known  enzymes.12 
Recent  research  interest  regarding  this  enzyme  is  not 
only  due  to  this  high  catalytic  efficiency  but  also  due  to 
the  broad  implications  of  AChE  inhibition  on  human 
health,  agrochemistry,  and  chemical  agents.  For  ex¬ 
ample,  Alzheimer's  disease  (AD)  is  associated  with  low 
in  vivo  levels  of  acetylcholine;  thus,  AChE  has  been 
targeted  in  many  drug  discovery  projects  aimed  at 
maintaining  ACh  availability  via  mild  or  reversible 
inhibitors  such  as  tacrine3  and  donepezil,4  etc.  While 
low-level  AChE  inhibition  is  useful  for  such  neurological 
treatments,  higher  levels  of  inhibition  can  be  detrimen¬ 
tal.  Organophosphorus  (OP)  compounds,  in  particular, 
irreversibly  deactivate  AChE  and  may  induce  failure  of 
cholinergic  synaptic  transmission,  deterioration  of  neu¬ 
romuscular  junctions,  flaccid  muscle  paralysis,  and 
central  nervous  system  seizures.5’6  Effective  drug  design 
thus  requires  great  care  in  balancing  the  level  of 
inhibitive  efficacy. 

The  availability  of  AChE  crystal  structures  for  various 
species  with  and  without  ligands  provides  a  solid  basis 
for  structure-based  design  of  novel  AChE  inhibitors.7 
There  aretwo  principle  binding  sites  in  theAChE.The 
catalytic  active  site  is  located  at  the  base  of  a  deep  gorge 
in  the  enzyme.  It  contains  a  catalytic  triad,  Ser203, 
Glu334,  and  His447  (huAChE  sequence  numbering,  to 
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be  used  throughout  unless  otherwise  specified),  and 
nearby  residues  (e.g.,  the  choline  binding  site:  Trp86) 
that  collectively  effect  the  ACh  catalysis  reactions.8 
AChE  also  has  a  peripheral  anionic  site  (PAS)  located 
near  the  enzyme  surface  at  the  mouth  of  the  active  site 
gorge.  The  residue Trp286  plays  a  very  important  role 
in  ligand  binding  in  the  PAS.  Ligand  binding  to  the  PAS 
affects  enzymatic  activity  through  a  combination  of 
steri  c  bl  ockade  of  I  i  gands  movi  ng  through  the  gorge  and 
allosteric  alteration  of  the  catalytic  triad  conformation 
and  efficiency.9 The  gorge  itself  is  a  narrow  hydrophobic 
channel  with  a  length  of  ~20A,  connecting  the  PAS  site 
to  the  active  site.10  An  acyl  binding  pocket  consists  of 
residues  Glyl22,  Trp236,  Phe295,  Phe297,  and  Phe338 
and  is  responsible  for  interacting  with  the  acetyl 
group.11  Early  inhibition  research  was  mainly  focused 
on  ligands  binding  in  the  active  site  (e.g.,  tacrine3, 
amiridine,12  etc.).  Recent  efforts  have  focused  on  finding 
novel  ligands  that  bind  to  both  sites  in  order  to  search 
for  more  potent  reversible  inhibitors  (e.g.,  TAK-147, 
E2020,  etc.),  selectivelyfavoringtheinhibition  of  AChE 
rather  than  the  related  butyryl cholinesterase  (BChE). 

Molecular  modeling  has  proven  increasingly  impor¬ 
tant  in  helping  to  design  novel  enzyme  inhibitors.  A 
substantial  amount  of  prior  AChE  inhibitor  research 
has  focused  on  using  ligand-based  design  methods  such 
as  CoM FA. 13-17  Given  an  accurate  receptor  structure, 
molecular  docking  can  also  be  very  useful  in  character¬ 
izing  ligand- receptor  binding  by  providing  predictions 
of  the  bound  conformation  for  the  ligand  and  a  scheme 
for  energetically  ranking  (i.e.,  scoring)  the  ligand- 
receptor  interaction.  Great  successes  have  been  achieved 
in  terms  of  conformational  predictions  via  flexible 
docking  programs  such  as  Dock,18  Gold,19  FlexX,20  etc. 
Such  conformational  predictions  are  very  important  to 
drug  design  because  (1)  the  binding  conformation  of 
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ligands  is  much  easier  to  validate  (i.e.,  through  com¬ 
parison  with  experimentally  observed  structures)  than 
is  the  binding  affinity  for  different  systems  and  (2)  the 
accurate  prediction  of  the  bound  conformation  is  a 
prerequisite  for  reliable  scoring.  Even  with  good  struc¬ 
tural  predictions,  however,  the  score  may  not  always 
agree  well  with  experimentally  determined  affinities 
mainly  because  experimental  conditions  include  impor¬ 
tant  dynamic  or  entropic  effects  that  are  difficult  to 
rigorously  represent  in  a  general  scoring  function.  To 
account  for  such  effects  empirically,  scoring  functions 
are  typically  trained  via  diverse  sets  of  previously 
characterized  ligand-receptor  interactions.  Unfortu¬ 
nately,  nofinite training  set  is  likely  to  provide  a  perfect 
representation  for  all  systems  of  interest  because  of  the 
varying  physicochemical  conditions  present  in  different 
receptors.  Indeed,  with  the  AChE  system,  studies  on 
steroidal  alkaloid  inhibition  show  no  correlation  be 
tween  the  calculated  binding  energy  and  experimentally 
determined  activity.21  Molecular  dynamics  simulations 
do  provide  a  natural  means  for  quantifying  both  the 
entropic  and  enthalpic  components  of  binding  affinity 
for  ligand-AChE  interactions;20  however,  extensive 
computational  demands  make  molecuar  dynamics  (M  D) 
simulations  prohibitively  time-consuming  for  analysis 
of  large  compound  collections.  In  such  cases,  the  best 
compromise  may  be  to  carry  out  simple  docking  studies 
but  to  explicitly  train  the  scoring  function  toreproduce 
behavior  in  the  system  of  interest.  Herein,  we  present 
a  new  scoring  function  that  is  based  on  the  ligand- 
receptor  interaction  field  and  is  trained  specifically  to 
reproduce  AChE  inhibition. 

Methods 

The  AChE  crystal  structures  used  herein  were  ob¬ 
tained  from  the  Protein  Databank  (PDB).  Our  main 
docking  and  training  activities  have  been  focused  on  a 
human  AChE  (huAChE)  structure  (code,  1B41)22  but 
rely  on  ligand  binding  information  from  a  Torpedo 
californica  AChE  (tcAChE)  structure  (code,  lEVE)that 
includes  a  cocrystallized  E2020  inhibitor.23  Given  the 
absence  of  a  firm  understanding  of  the  persistence  and 
roles  of  individual  solvent  molecules  in  and  around  the 
AChE  binding  sites,  all  waters  were  removed  from  the 
structures.  To  ascertain  the  orientation  of  the  ligand 
E2020  relative  to  the  huAChE  structure,  the  huAChE 
structure  was  aligned  to  the  tcAChE  in  SYBYL24  by 
achieving  a  maximal  overlap  of  Ca  atoms  for  corre¬ 
sponding  huAChE/tcAChE  residues  within  the  recep¬ 
tor  region.  The  resulting  root-mean  square  deviation 
(rmsd)  between  thetwo  aligned  huAChE/tcAChE  struc¬ 
tures  is  only  0.85  A  for  the  set  of  all  backbone  Ca  atoms 
within  the  full  enzyme  subunit,  suggesting  good  overall 
alignment  and  substantial  structural  similarity.  Hy¬ 
drogen  atoms  were  added  (via  SYBYL)  to  the  resulting 
huAChE-E2020  complex.  The  positions  of  these  new 
protons  were  then  optimized  in  MOE25  via  molecular 
mechanics  using  the  MMFF94s  force  field26  (all  heavy 
atoms  fixed)  to  avoid  bad  interatomic  contacts.  The 
position  of  E  2020  was  then  optimized  (all  receptor  atoms 
fixed)  to  determine  a  plausible  stable  conformational 
structure  for  the  ligand  in  the  receptor  environment. 
In  both  of  the  above  simulations,  MMFF94s  charges 
were  used  to  account  for  relevant  electrostatics.  The 


steepest  descent  minimization  algorithm  was  used  for 
the  first  100  steps  (unless  an  rms  gradient  of  less  than 
100  kcal/(mol-A)  was  first  achieved),  followed  by  200 
steps  of  conjugate  gradient  (unless  an  rms  gradient  of 
less  than  1  kcal/(mol*A)  was  attained),  and  finally 
completed  by  1000  steps  of  truncated  Newton  (or  an  rms 
gradient  of  less  than  0.01  kcal/(mol-A)).  The  resulting 
E2020  structure  was  then  extracted  for  subsequent 
docking  calculations. 

Sixty-nine  compounds  with  ICsodata  measured  with 
human  AChE  assay27-30  were  selected  for  training  and 
testing  the  scoring  function.  The  activity  among  these 
compounds  ranges  from  0.33  to  30  000  nM  (T ables  1  and 
2).  The  active  site  for  the  huAChE  docking  calculations 
was  constructed  from  the  crystal  structure  by  retai  ni  ng 
all  residues  within  a  radius  of  12  A  relative  to  E2020 
(but  discarding  the  original  ligand  itself).  Docking 
calculations  were  carried  out  with  the  Gold  program.19 
A  genetic  algorithm  was  used  in  searching  the  binding 
conformation  of  flexible  ligands,  using  the  default 
parameters  in  GOLD.  A  maximum  of  20  poses  were 
computed  for  each  compound.  Those  docked  conforma¬ 
tions  were  saved  in  SDF  format  and  then  imported  into 
SYBYL  for  scoring  calculations  according  to  the  FlexX 
and  CSCORE  modules.  The  scoring  methods  available 
included  empirical  methods  such  asChemScore,31  FlexX 
score,20  and  G  Score19  and  knowledge-based  methods 
such  as  PMF  score32  and  DrugScore.33  Multilinear 
regressi on  (M  L  R)  was  used  to  obtai n  a  consensus  score 
from  these  methods.  One  conformation  was  selected  for 
each  compound  to  give  a  good  compromise  between  the 
best  consensus  score  and  those  with  the  closest  align¬ 
ment  to  the  original  E 2020  ligand.  Specifically,  the  pose 
for  the  scoring  of  the  activity  was  selected  on  the  basis 
of  having  the  highest  consensus  score  (first  criterion) 
and  ChemScore  (tie-breaking  criterion),  with  thefurther 
stipulation  that  thefollowing  knowledge-based  criteria 
(as  determined  by  visual  inspection)  must  be  obeyed 
whenever  possible:  (1)  good  TV  TV  overlap  with  residue 
Trp86,  as  has  been  found  to  be  critical  for  E2020 
binding;23  (2)  good  TV  TV  overlap  with  residue  Trp286, 
as  has  also  been  found  to  be  very  important  for  E2020 
complexation23 

The  chosen  conformations  were  used  to  fit  an  interac¬ 
tion  field  whose  form,  basically  a  variant  of  the  com¬ 
parative  binding  energy  (COM  Bl  N  E )  method,34-36  is  as 
fol  I  ows: 

pi  Cm  =  £c,e«  +  ■£<!&»  (l) 

i  J 

where  q  and  dj  are  fitted  coefficients,  Efle  and  £jdw  are 
electrostatic  and  van  der  Waals  (VDW)  interactions 
arising  between  atoms  in  the  / th  and y'th  residues  and 
the  ligand.  In  this  expression,  all  receptor  residues 
within  10  A  of  the  position  of  the  original  E 2020  ligand 
(a  total  of  92  residues)  wereinduded  in  thesummation 
over  /,  and  E?e  and  E/dw  were  calculated  via  an  SVL 
script  written  for  the  MOE  system.  The  statistic  analy¬ 
sis  was  performed  in  Simca-P37  with  partial  least 
squares  regression  (PLS).  Fifty-three  of  the  full  69 
compounds  were  selected  as  our  training  set,  and  the 
other  16  compounds  were  used  as  a  test  set  for  validat¬ 
ing  the  predictive  power  of  the  new  scoring  function. 
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Results  and  Discussion 

Our  scoring  model  built  via  PLS  regression  over 
i nteracti ons  wi thi n  the  53-molecul e trai ni ng  set  appears 
to  be  of  reasonable  quality,  with  a  correlation  coefficient 
of  R2  =  0.89  and  a  leave-one-out  cross-validation  cor¬ 
relation  of  Q2  =  0.72.  In  using  the  scoring  function  to 
evaluate  the  activity  of  the  16-molecule  testing  set,  we 
achieved  good  predictivity:  a  correlation  of  R2  =  0.69 
(Figure  1)  between  the  calculated  results  and  experi¬ 
mental  values. 


To  compare  the  precision  and  extensibility  of  our 
scoring  function,  we  contrasted  the  above  results  with 
predictions  made  using  several  commercially  available 
scoring  methods,  including  ChemScore,31  FlexX  score,20 
DrugScore,33  G  Score,19  and  PM  F  score.32  The  correla¬ 
tion  between  the  experiment  and  any  single  scoring 
method  is  poor.  The  PMF  score  showed  the  best  cor¬ 
relation  but  was  still  poor  (R2  =  0.13  for  the  training 
set).  All  of  the  other  representations  gave  even  worse 
correlations:  ChemScore  =  0.07,  FlexX  =  0.05,  Drug- 
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Correlation  for  the  training  set 


Correlation  for  the  testing  set 
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Figure  1.  Correlation  of  the  calculated  activity  (pi C5o)  with 
experiment:  (top)  training  set  (R2  =  0.89);  (bottom)  testing 
set  (R2  =  0.69). 

Score  =  0.01,  G  Score  =  0.004.  Since  those  general 
scoring  methods  were  not  trained  in  this  AChE  inhibi¬ 
tion  set,  it  is  not  completely  fair  to  compare  them 
directly  with  our  specially  tailored  energy  decomposition 
method.  Therefore,  we  built  a  consensus  score  by 
training  (over  the  53-molecule  set)  a  weighted  sum  over 
the  above  five  commercially  avail  able  scoring  functions 
as  follows: 

pi  C50  =  0.05682  Chemscore  -  0.00499  Drugscore  - 
0.03582  Flexscore-  0.01232  Gscore  + 

0.01847  PM  F score  +  7.62820  (2) 

where  Chemscore,  Drugscore,  Flexscore,  Gscore,  and 
PMFscore  refer  to  computed  affinities  from  the  Chem¬ 
score,  DrugScore,  FlexX  score,  G  Score,  and  PMF  score 
methods,  respectively.  Although  an  improvement  was 
found  for  this  consensus  score  within  the  training  set 
itself  (R2  =  0.26),  its  predictive  potency  is  poor,  judging 
by  no  evident  correlation  within  the  16-moleculetesting 
set. 

To  help  verify  the  physical  sensibility  of  our  model, 
we  have  mapped  out  the  residues  that  contribute 
significantly  to  the  scoring  function.  The  coefficients  of 
the20  most  important  residues  in  terms  of  electrostatic 
and  VDW  contributions  are  shown  in  Figure  2.  In  those 
residues,  T rp86, 1  Ie451,  Gly448,  Tyr449,  and  Ser229 are 
the  most  important  residues  in  the  active  site  for  VDW 
interactions.  Trp86  functions  by  forming  , 7t  JZ  interac¬ 
tion  with  ligand  aryl  groups  (when  available),  whilethe 
other  residues  define  the  shape  of  the  gorge  base, 
servi  ng  to  discri  mi  nate  accordi ng  to  I igand  shape.  I  n  the 
upper  gorge  area  and  the  acyl  binding  pocket,  residues 
Tyrl24,  Phe295,  Phe338,  and  Phe297  are  responsible 
for  providing  hydrophobic  contacts.  The  ring  of  Tyr72 


Residue  importance  for  electrostatic 
interactions 


Residue  importance  for  VDW  interactions 


Residues 


Figure  2.  (Top)  Coefficients  of  the  20  most  important  residues 
for  electrostatic  interactions.  (Bottom)  Coefficients  of  the  20 
most  important  residues  for  VDW  interactions. 

is  almost  perpendicular  totheTrp286  ring  and  forms  a 
blocking  wall  to  prevent  the  ligand  ring  from  moving 
away  from  the  position  where  it  forms  a  JZ~JZ  interaction 
with  theTrp286  ring.  Phe295,  Phe297,  Val365,  and 
Glu292  form  another  wall  on  the  other  side  of  the  gorge, 
stretching  from  the  acyl  pocket  toward  the  PAS. 

Residues  Tyr449,  Glu450,  Ile451,  Alal27,  Serl28, 
Tyrl33,  Ilell8  near  Trp86,  and  the  "oxyanion  hole" 
residues  Glyl21  and  Glyl22  are  important  in  providing 
electrostatic  interactions  in  the  active  site.  Tyr337, 
Asp74,  Thr83,  and  Asn87  are  the  primary  electrostatic 
contributors  in  the  gorge  area.  Gly342,  Leu 76,  Glu285, 
Trp286,  His287,  and  Gln291  are  probably  helpful  in 
enhancing  the  activity  of  ligands  with  polar  groups 
oriented  in  this  area,  as  is  evidenced  by  reports  that  an 
AChE  inhibitor  tethering  in  the  position  of  His287  can 
affect  the  binding  affinity  as  much  as  14-fold.38  Site- 
directed  mutagenesis  in  huAChE  indicated  that  Asp74, 
Tyr337,  Phe338,  Phe295,  Phe297,  Tyrl33,  and  Glu450 
can  affect  the  affinity  although  it  has  been  difficult  to 
determine  experimentally  whether  these  residues  con¬ 
tribute  mainly  electrostatic  or  VDW  interactions.39’40 

The  docked  ligand  structures  generally  support  the 
above  analysis  regarding  the  identity  of  principle  resi¬ 
dues.  In  our  current  docking  calculations,  molecules  2-7 
all  share  similar  conformations  in  the  PAS.  The  modi¬ 
fied  groups  in  those  molecules  are  actually  exposed  to 
thesolvent  and  do  not  contribute  directly  to  the  I  igand- 
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Figure  3.  (T op)  Structure  of  E 2020  bi  ndi  ng  to  tcAChE .  E 2020 
is  rendered  in  cyan.  (Bottom)  Structure  of  the  E 2020  binding 
in  huAChE.  E2020  is  rendered  in  cyan. 


receptor  interaction.  This  is  the  same  as  observed  in 
previous  studies.28  However,  in  the  current  work,  mol¬ 
ecule  8  takes  a  slightly  different  conformation  with  its 
morpholi no  group  situating  in  the  half-buried  pocket  by 
Trp286,  His287,  Ser293,  Glu292,  and  Leu289.  The 
morpholi  no  oxygen  has  a  distance  of  3.65  A  from  the 
backbone  N  of  Glu292,  and  the  morpholino  nitrogen  is 
3.78  A  from  the  backbone O  of  Ser 293.  Such  interactions 
could  slightly  pull  the  benzisoxazole  ring  away  from 
Trp286  ring,  leaving  only  a  partial  jv—jv  interaction.  This 
particular  conformation  leads  to  an  affinity  increased 
by  more  than  10-fold  to  0.8  nM  relative  to  molecules 
2-7.  This  particular  region  has  been  confirmed  in  an 
X-ray  structure  to  be  very  important  for  the  inhibitors 
binding  to  the  PAS.9 

As  a  final  point  of  validation,  we  compared  the 
calculated  conformation  for  E2020  within  the  huAChE 
crystal  structure  relative  to  its  original  cocrystallized 
conformation  in  tcAChE.  Our  calculated  structure  sug¬ 
gests  that  E2020  has  similar  but  not  identical  binding 
modes  in  tcAChE  and  huAChE  (Figure3).  In  theactive 
site,  its  benzyl  ring  forms  a  TZ  71  interaction  with  the 
indole  ring  of  Trp86  in  huAChE  and  Trp84  in  tcAChE. 
In  the  PAS,  the  indanone  ring  of  E2020  forms  a  7Z—7T 
interaction  with  the  indole  ring  of  Trp286  in  huAChE 
and  Trp279  in  tcAChE.  The  charged  nitrogen  of  the 
E2020  pi  peri  dine  ring  undergoes  a  cation-jr  interaction 
with  the  phenyl  ring  of  Phe330  in  tcAChE.  The  corre¬ 
sponding  residue  Tyr337  in  huAChE  does  not  form  a 
similar  cation-jr  interaction  due  to  the  steric  I  imitations 


in  this  area  associated  with  an  inauspicious  orientation 
of  the  Tyr337  ring.  As  a  result,  the  nitrogen  of  the 
piperidine  ring  of  the  E2020  instead  interacts  with  the 
hydroxyl  groups  of  the  Tyr337,  Tyrl24,  and  Ser  125 
within  distances  of  3.41,  3.12,  and  4.18  A  relative  to  the 
oxygen  atom  of  the  respective  hydroxyl  groups.  T  o  probe 
the  role  of  Tyr337,  we  examined  the  potential  energy 
curve  for  Tyr 337  side  chain  rotation  relative  to  theother 
resi  dues  (energy  eval  uati  on  accordi  ng  to  the  M  M  F  F  94s 
force  field  with  appropriate  charges)  but  found  only  one 
minimum  in  the  potential  curve.  Closer  analysis  reveals 
that  that  theTyr337  ring  is  trapped  in  a  local  pocket 
formed  by  Phe338,  Tyr 341,  Trp439,  Trp449,  and  His447. 

I  n  previous  structural  studies  and  molecular  dynamics 
simulations,  it  has  been  found  that  Phe330  in  tcAChE 
can  adopt  a  wide  range  of  conformations  in  the  complex 
structure  and  may  function  as  a  swinging  gate1523'39 
that  structural  ly  couples  the  anionic  subsite  of  the  active 
site  and  the  PAS.  It  is  natural  to  expect  similar  behavior 
of  Tyr 337  in  the  huAChE  compared  with  the  analogous 
Phe330  in  tcAChE.  Such  a  gate  swing  movement  of  the 
Tyr 337  ring  would  certainly  induce  a  shift  in  attached 
and  neighboring  residues;  thus,  one  function  of  this 
PAS-active  site  coupling  could  be  to  effect  proper 
residue  alignment  within  the  anionic  subsite.  The  fact 
that  the  huAChE  crystal  structure  used  in  this  study 
to  train  our  scoring  function  did  not  have  a  cocrystal¬ 
lized  active  site  inhibitor  (having  only  a  PAS-bound 
fasciculin  molecule22),  and  thus  did  not  reflect  such  a 
conformational  shift  in  theTyr337,  may  constitute  a 
subtleflaw  in  our  model.  Given  our  model's  fairly  strong 
predictive  capacity,  we  expect  the  flaw  to  be  of  only 
minor  consequence,  however. 

The  accord  between  results  derived  from  our  scoring 
method  and  those  of  X-ray  structures  and  mutagenesis 
observations  indicates  the  effectiveness  of  the  current 
analysis  and  the  scoring  function's  predictive  power. 
Since  this  method  requires  a  set  of  compounds  with 
known  activity  data  to  fit  the  score  function,  it  is  most 
applicable  to  the  task  of  lead  optimization  as  opposed 
to  de  novo  discovery.  I  n  contrast  to  the  comparative 
molecular  field  analysis  (CoMFA)41  method,  which  uses 
probe  atoms  such  as  nitrogen,  carbon,  etc.  to  calculate 
possible  interaction  fields  between  ligands  and  a  puta¬ 
tive  receptor  (and  does  not  require  explicit  atomic-level 
representation  of  the  receptor  structure),  our  method 
uses  a  "real"  interaction  field  between  the  ligands  and 
receptor,  thus  requiring  advance  knowledge  of  the 
receptor  structure.  The  benefit  of  using  the  current 
method  is  that  the  important  residues  can  be  directly 
identified  and  verified  by  site-directed  mutagenesis, 
whereas  CoMFA  merely  generates  a  hypothetical  map 
of  possible  favorable  and  unfavorable  interaction  regions 
based  on  statistically  postulated  correlations  between 
activity  and  orientations  of  specific  functional  groups 
within  the  ligand  training  set.  Such  a  hypothetical  map 
may  correspond  to  the  real  interaction  field  if  the 
underlying  postulated  statistical  correlation  is  valid  but 
is  unlikely  to  be  as  accurate  a  representation  as  a  real 
interaction  field.  In  cases  where  crystal  structures  of 
the  relevant  receptor  are  available,  CoMFA-generated 
hypotheses  may  be  adjusted  in  order  to  directly  cor¬ 
respond  to  the  real  field.  I  n  this  vein,  there  have  been 
successful  studies  that  used  receptor-based  docking 
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methods  to  al  i  gn  the  I  i  gands  and  then  used  CoM  F  A  on 
the  resulting  alignments  to  ascertain  and  analyze  the 
resulting  grid  for  drug  design  purposes.13-14'16-17  This 
method  is  similar  in  principle  to  the  COMBI NE-like 
technique  that  we  have  applied  herein  to  this  huAChE 
system  except  that  COMBINE  methods  are  easier  to 
compute  and  to  understand  because  they  only  require 
treatment  of  ligand-residue  interactions  (hundreds  of 
terms  or  less)  rather  than  the  thousands  of  grid  points 
in  CoM  FA  models.  Another  benefit  relative  to  CoM  FA 
models  is  that  the  relative  simplicity  of  our  scoring 
function  also  allows  it  to  be  used  directly  in  any 
subsequent  docking  calculations,  thus  providing  a  means 
to  accurately  score  novel  inhibitor  candidates  and 
predict  their  bound  conformations  rather  than  having 
to  do  so  in  a  stepwise  manner. 

Conclusions 

A  new  AChE-specific  scoring  function  has  been  de¬ 
veloped  herein  and  used  in  predicting  binding  affinity 
of  AChE  inhibitors.  The  method  is  based  on  a  COM¬ 
BI  N  E  -type  approach,  i  ncorporati  ng  the  el  ectrostati c  and 
VDW  interaction  fields  between  the  ligands  and  recep¬ 
tor.  A  53-compound  training  set  was  used  to  construct 
the  scoring  function,  and  a  further  16  compounds  were 
used  to  test  the  resulting  model.  Strong  statistical 
correlations  were  found  between  predicted  and  observed 
affinities  for  both  thetraining  and  testing  set.  Analysis 
of  the  scoring  function  has  permitted  identification  of 
those  receptor  residues  making  the  most  important 
contributions  to  ligand  binding.  These  analyses  are 
consistent  with  the  X-ray  structure  and  mutagenesis 
studies.  A  comparison  with  other  scoring  methods  and 
consensus  scoring  indicated  the  high  effectiveness  and 
predictability  of  this  method. 
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Acute  toxicity  of  organophosphorus  (OP)  compounds  results  mainly  from  irreversible  acetylcholinest¬ 
erase  (AChE)  inhibition;  however  OP  toxicity  frequently  hinges  on  prior  biotransformations  that  produce 
toxic  metabolites.  To  account  for  both  precursor  metabolic  effects  and  primary  AChE  inhibition,  we 
included  absorption,  distribution,  metabolism,  excretion  (ADME)  effects,  ligand  binding,  and  reactive 
AChE  phosphorylation  and  aging  in  a  detailed  but  computationally  expedient  phenomenological  toxicity 
model.  Ligand  negative  accessible  surface  area  (NASA)  was  used  as  a  generic  ADME  descriptor,  while 
relevant  metabolic,  phosphorylation,  and  aging  reactions  were  assessed  via  quantum  chemical  enthalpy 
calculations,  and  the  binding  affinity  of  the  Michaelis  complex  was  quantified  via  Comparative  Molecular 
Field  Analysis  (CoMFA).  The  resulting  model  correlates  very  well  ( R 2  =  0.90)  with  experimental  acute 
toxicity  measurements  and  provides  useful  mechanistic  insight  into  the  underlying  toxicity.  Model 
predictivity  was  validated  by  leave-one-out  cross-validation  (Q2  =  0.82).  The  Michaelis  binding  affinity 
descriptor  has  the  largest  weight  in  our  model,  but  subsequent  covalent  inhibition  and  prior  ADME  effects 
also  exhibit  significant  effects. 


Introduction 

Organophosphorus  (OP)  compounds  are  best  known  for  their 
highly  toxic  effect  as  agricultural  chemicals  and  chemical 
warfare  agents.  The  mechanism  of  the  acute  OP  toxicity  has 
been  a  subject  of  substantial  research  interest  for  several  decades 
in  the  search  for  therapies  to  counter  OP  poisoning  (e.g.,  oximes 
(7)),  to  devise  more  effective  pesticides,  and  to  uncover  low- 
toxicity  OP  species  for  potential  pharmacological  applications 
(2).  The  main  cause  of  OP  toxicity  is  irreversible  inhibition  of 
the  acetylcholinesterase  (AChE)  enzyme  (2).  AChE  hydrolyzes 
the  neurotransmitter  acetylcholine  (ACh)  at  cholinergic  syn¬ 
apses.  The  inactivation  of  the  AChE  by  OP  compounds  thus 
causes  acetylcholine  buildup,  which  can  in  turn  induce  failure 
of  cholinergic  synaptic  transmission,  leading  to  deterioration 
of  neuromuscular  junctions,  flaccid  muscle  paralysis,  and  central 
nervous  system  seizures  ( 3 ,  4). 

Numerous  studies  via  a  number  of  different  techniques  (e.g., 
crystallography  (5),  mutagenesis  (<5),  modeling  (7),  etc.)  concur 
that  OP  compounds  mainly  interact  with  the  AChE  catalytic 
triad  (residues  Ser203,  His447,  and  Glu334  according  to  the 
mouse  AChE  sequence)  in  the  active  site,  which  is  connected 
by  a  deep  narrow  gorge  to  the  peripheral  binding  site  at  the 
entrance  of  the  enzyme.  The  inhibition  mechanism  involves 
phosphorylation  whereby  a  covalent  bond  between  the  central 
phosphorus  atom  of  the  ligand  and  the  side  chain  oxygen  of 
Ser203  is  formed.  For  some  OP  compounds,  the  resulting 
conjugate  can  further  react  via  a  process  known  as  “aging”  that 
usually  involves  the  loss  of  an  alkyl  group  from  the  phosphyl 
alkoxy  substitutent.  The  aging  process  prevents  AChE  reactiva¬ 
tion  by  conventional  antidotes  such  as  oxime  therapy. 
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Many  OP  compounds  do  not  interact  with  AChE  in  vivo  in 
their  native  form.  The  metabolite  products,  usually  oxon 
analogues,  are  frequently  the  actual  toxicants  ( 8 ).  The  usual 
metabolism  reactions  include  oxidative  desulfuration,  N-dealky- 
lation,  O-dealkylation,  O-dearylation,  thioether  oxidation,  among 
others  (9).  It  is  thus  not  surprising  that  in  vitro  inhibition 
measurements  rarely  provide  a  simple  linear  correlation  with 
the  acute  toxicity.  To  better  understand  the  mechanism  of  OP 
activity,  it  is  necessary  to  augment  information  from  direct 
AChE  inhibition  studies  with  an  account  of  absorption,  distribu¬ 
tion,  metabolism,  and  excretion  (ADME)  processes.  Given  such 
complex  underlying  biological  processes,  a  full  investigation 
and  understanding  of  the  uptake  and  toxicology  of  even  a  single 
OP  compound  can  make  for  a  daunting  experimental  challenge. 

Molecular  modeling  methods  have  been  very  useful  in  helping 
to  understand  OP  toxicity  because  of  the  complexity  of  the 
AChE  system  and  the  arduous,  expensive,  and  sometimes  dan¬ 
gerous  nature  of  the  corresponding  experiments.  Previous  model¬ 
ing  studies  have  largely  focused  on  the  direct  inhibition  of 
AChE,  mainly  by  probing  the  interaction  of  ligands  with  the 
AChE  active  site  (10).  Those  methods  include  quantum  chem¬ 
istry  (7),  molecular  dynamics  (7),  docking  (77),  and  interaction 
field  analysis  of  noncovalent  ligands  (72).  In  the  present  paper, 
we  will  focus  on  a  detailed  but  computationally  efficient  scheme 
for  integrating  ADME  properties,  AChE  binding  affinity, 
phosphorylation,  and  the  aging  process  into  a  single,  unified, 
phenomenological  model  capable  of  reproducing  the  primary 
trends  in  OP  toxicity  in  a  mechanistically  intuitive  fashion. 

Materials  and  Methods 

Inhibitors.  Data  corresponding  to  the  38  OP  compounds 
considered  for  this  study  are  listed  in  Table  1.  The  LD50  values 
correspond  to  studies  on  male  rats  exposed  orally  and  are  given  in 
mg/kg  body  weight.  Among  those  compounds,  Dichlorvos,  Dicro- 
tophos,  Naled,  Tetrachlorvinphos,  Ethoprop,  Oxydemeton-methyl, 
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Table  1.  Organophosphorus  Compounds  Considered  for  Analysis" 


parent 

metabolite 

direct 

acting 

metabolite 

references 

LDso" 

(mg/kg) 

pLD5oc 

Dichlorvos 

yes 

80 

4.097 

Dicrotophos 

yes 

21 

4.678 

Naled 

yes 

250 

3.602 

T  etrachlorvinpho  s 

yes 

4000 

2.398 

Trichlorfon 

Dichlorvos 

no 

31 

650 

3.187 

Ethoprop 

yes 

34 

4.469 

Azinphos  Methyl 

oxon 

no 

32 

13 

4.886 

Bensulide 

? 

no 

770 

3.113 

Dimethoate 

oxon 

no 

33 

215 

3.668 

Disulfoton 

? 

no 

2 

5.699 

Ethion 

oxon 

no 

34 

13 

4.886 

Malathion 

oxon 

no 

35 

1375 

2.862 

Methidathion 

oxon 

no 

36 

31 

4.509 

Phorate 

oxon 

no 

37 

2 

5.699 

Phosmet 

oxon 

no 

38 

147 

3.833 

Sulfopros 

oxon 

no 

39 

65 

4.187 

Temephos 

? 

no 

8600 

2.066 

Terbufos 

oxon 

no 

41 

2 

5.699 

Fonofos 

oxon 

no 

41 

3 

5.523 

Oxydemeton-methyl 

yes 

42 

47 

4.328 

Profenofos 

yes 

43 

358 

3.446 

Chlorethoxyfos 

? 

no 

5 

5.301 

Chlorpyrifos 

oxon 

no 

44 

155 

3.810 

Coumaphos 

oxon 

no 

45 

41 

4.387 

Diazinon 

oxon 

no 

46 

250 

3.602 

Fenitrothion 

oxon 

no 

47 

740 

3.131 

Fenthion 

oxon 

no 

48 

215 

3.668 

Methyl  Chlorpyrifos 

oxon 

no 

49 

1500 

2.824 

Methyl  parathion 

oxon 

no 

50 

14 

4.854 

Parathion 

oxon 

no 

51 

13 

4.886 

Pirimiphos  methyl 

oxon 

no 

52 

1450 

2.839 

Sulfotep 

? 

no 

5 

5.301 

Tebupirimphos 

? 

no 

2 

5.699 

Fenamiphos 

yes 

8 

5.097 

Acephate 

? 

no 

700 

3.155 

Methamidophos 

yes 

31 

4.509 

Isofenfos 

oxon 

no 

53 

28 

4.553 

Propetemphos 

? 

no 

119 

3.924 

a  The  compounds  with  question  mark  in  the  metabolite  were  not  used  in 
the  current  work.  b  The  data  are  from  ref  8.  c  The  pLDso  was  calculated  as 
6  -  log(LD50). 


Profenofos,  Fenamiphos,  and  Methamidophos  all  have  direct 
interaction  with  AChE,  while  the  other  29  undergo  biotransforma- 
tion  before  reaching  the  AChE  active  site.  Unfortunately,  to  the 
best  of  our  knowledge,  there  are  no  conclusive  metabolic  data 
available  in  the  literature  for  Bensulide,  Disulfoton,  Temephos, 
Chlorethoxyfos,  Sulfotep,  Tebupirimphos,  Acephate,  and  Pro- 
peteamphos.  Based  on  data  availability,  we  thus  chose  a  total  of 
30  OP  compounds  (9  direct-acting  and  21  indirect- acting  species 
for  which  definitive  metabolic  data  exists)  for  further  analysis. 

Receptor  Model.  The  following  AChE  crystal  structures  were 
obtained  from  protein  data  bank  (13)  and  used  as  structural 
references  in  this  study:  (1)  Torpedo  calif ornica  AChE  with  nerve 
agent  VX  after  phosphorylation  (1VXR)  (14),  (2)  Torpedo  Cali¬ 
fornia  AChE  with  nerve  agent  VX  after  aging  (1VXO)  (14),  and 
(3)  mouse  AChE  complexed  with  Gallamine  (1N5M)  (15).  The 
latter  structure,  with  Gallamine  removed,  was  used  as  a  template 
for  assembling  our  ligand— receptor  complex  models.  Since  rat 
AChE  is  highly  similar  to  mouse  AChE  (only  11  amino  resides 
difference  between  the  two  (16),  with  none  of  those  11  playing 
any  known  role  in  the  enzymatic  activity  or  inhibition),  we  expect 
no  appreciable  difference  in  result  to  arise  from  using  the  mouse 
rather  than  the  rat  structure.  All  water  and  heteroatoms  were 
removed  from  the  structure  before  adding  hydrogen  atoms.  Protein 
electrostatics  were  modeled  via  AMBER94  charges  (17),  and 
MMFF94  charges  (18)  were  employed  for  ligands.  The  1VXR  and 
1 VXO  structures  served  as  guides  for  positioning  inhibitors  within 
the  AChE  structure.  Both  1VXR  and  1YXO  were  superimposed 
with  1N5M  to  determine  (by  analogy  to  the  cocrystallized  VX) 
positions  of  the  ligands  in  pro-aging  and  post-aging  by  methods 
that  have  been  described  in  a  previous  study  on  noncovalent 
inhibitors  (12). 


Descriptors.  Given  a  focus  on  developing  a  mechanistically 
insightful  phenomenological  model,  descriptors  were  chosen  ac¬ 
cording  to  a  requirement  of  having  a  natural  physicochemical 
relationship  with  the  various  physiological  effects  under  consider¬ 
ation  (ADME  processes,  noncovalent  complex  formation,  and 
covalent  reactivity).  For  reasons  elaborated  on  in  the  coming 
sections,  ADME  effects  were  modeled  via  ligand  negative  accessible 
surface  area  (NASA;  as  computed  via  the  MOE  suite  of  software 

(19)  using  default  settings  and  MMFF94  charges  (18))  plus 
thermochemical  characterization  of  the  associated  metabolic  reac¬ 
tion,  noncovalent  binding  was  approximated  via  a  CoMFA  model, 
and  trends  among  covalent  enzymatic  and  inhibitive  processes  were 
modeled  via  their  thermochemistry.  All  thermochemical  quantifica¬ 
tion  was  done  via  quantum  mechanically  computed  molecular 
energies.  To  strike  a  balance  between  computational  efficiency  and 
reliable  representation  of  energetic  trends  among  our  ligands, 
quantum  chemical  calculations  were  carried  out  at  the  AMI  level 

(20)  using  MOPAC7  (27).  The  MMFF94  force  field  (18)  was  used 
in  the  various  nonquantum  chemical  molecular  modeling  calcula¬ 
tions.  Conformation  searches  were  performed  using  a  stochastic 
method  similar  to  the  RIPS  method  (22).  Ligand  alignment  for  the 
CoMFA  model  was  determined  by  a  pharmacophore-driven  flexible 
alignment  method  within  MOE  (19).  In  the  interaction  field 
calculations,  a  monoanionic  oxygen  anion  was  used  as  probe  atom 
to  calculate  the  energy  on  a  grid  at  intervals  of  1  A.  Both 
electrostatic  and  van  der  Waals  interaction  energies  were  saved  at 
each  grid  point  for  further  partial  least-squares  (PLS)  analysis  via 
the  Simca-P  software  (23). 

Results  and  Discussion 

Metabolic  Effects.  OP  compounds  may  be  biotransformed 
by  many  of  the  different  cytochrome  P450  isozymes.  For 
example,  CYP2C11,  CYP3A4,  and  CYP2B1/2  were  recently 
identified  as  playing  a  role  in  rat  metabolism  of  diazinon  (9). 
At  this  current  stage,  there  is  little  information  on  which  specific 
isozyme  accomplished  the  metabolism  of  each  specific  OP 
compounds.  In  general,  the  reaction  rate  of  the  biotransformation 
should  be  determined  largely  by  the  reaction  barrier  of  the 
corresponding  P450  enzymatic  reaction.  Theoretically,  these 
barriers  could  be  accurately  quantified  via  high-level  quantum 
mechanical  transition  state  calculations,  however  such  a  com¬ 
putationally  demanding  undertaking  is  not  realistic  given  the 
number  of  OP  compounds  to  be  studied  and  the  inherent 
structural  complexity  of  enzymatic  reactions,  not  to  mention 
the  frequent  lack  of  information  regarding  the  specific  P450 
isozyme  responsible  for  metabolizing  a  given  compound. 
However,  it  is  substantially  easier  to  identify  and  characterize 
specific  metabolites  rather  than  the  intervening  transition  states. 
As  a  simple  approximation,  we  thus  consider  the  underlying 
thermochemistry,  energetically  characterizing  both  the  initial 
and  final  states  for  the  biotransformation,  since  those  two  states 
do  carry  important  information  about  the  viability  and  propensity 
of  a  given  biotransformation.  All  relevant  initial  and  final  state 
compounds  are  listed  in  Table  1,  under  the  “parent”  and 
“metabolite”  columns,  respectively.  All  of  the  toxic  OP  com¬ 
pounds  listed  in  Table  1  were  identified  from  a  literature  search 
on  rat  OP  toxicology.  Most  of  the  metabolites  take  the  form  of 
the  corresponding  OP— oxon.  Trichlorfon  has  been  biotrans¬ 
formed  as  Dichlorvos,  which  itself  can  have  a  direct  interaction 
with  AChE.  Interestingly,  Trichlorfon  has  a  LD50  of  8  times 
greater  than  its  metabolite  Dichlorvos,  a  difference  likely 
reflecting  the  dilutional  effects  of  the  biotransformation. 

Molecular  Surface  Area.  The  surface  area  has  been  found 
to  play  important  roles  in  various  ADME  processes  for  central 
nervous  system  drugs,  including  effects  on  membrane  transport 
(24),  Caco-2  permeability  (25),  and  blood— brain  barrier  transit 
(26),  among  others.  Since  the  negative  and  positive  surfaces 
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Figure  1.  Correlation  of  the  negative  polar  surface  area  between  parent 
OP  and  its  metabolite. 
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Figure  2.  Energy  and  surface  area  of  the  3000  conformations  for 
ethion.  The  filled  circles  indicate  the  100  conformations  with  lowest 
energy. 

play  different  roles  in  the  ADME  properties,  it  is  necessary  to 
deal  with  them  separately.  In  OP  compounds,  one  of  the  most 
obvious  distinguishing  structural  characteristics  is  the  distribu¬ 
tion  of  negative  charge  in  the  shell  around  the  phosphate, 
phosphonate,  or  phosphinate  group.  Consequently,  we  have 
selected  the  negative  accessible  surface  area  (NASA)  as  a 
primary  descriptor  to  account  for  the  ADME  properties  of  OP 
compounds.  Since  both  the  parent  OP  and  its  metabolite  undergo 
some  distribution  and  transportation,  these  properties  should  be 
considered  for  estimating  ADME  properties.  We  compared  the 
NASA  between  parent  OP  compounds  and  their  metabolites 
(lowest  energy  conformations  in  each  case)  and  report  the  results 
in  Figure  1.  A  strong  correlation  (. R 2  =  0.98)  is  observed 
between  parent  and  metabolite  NASA,  which  is  not  surprising 
because  the  OP— oxon  form  of  most  of  the  metabolites  is  highly 
similar  to  the  parent  structure.  In  this  regard,  NASA  values 
derived  from  either  the  parent  OP  or  its  metabolite  are  equally 
suitable  for  estimating  ADME  effects. 

Since  NASA  values  depend  on  the  molecular  conformation, 
it  is  necessary  to  take  into  account  any  conformation  change 
that  may  result  from  constituent  steps  in  the  ADME  process. 
Conformational  flexibility  is  related  to  the  number  of  the 
rotatable  bonds  within  a  molecule.  In  the  30  OP  compounds, 
the  number  of  rotatable  bonds  ranges  from  5  for  methamidophos 
to  16  for  ethion.  Across  a  spread  of  3000  distinct  conformations 
of  ethion  (Figure  2),  the  NASA  varies  by  at  most  20%;  however, 
in  the  100  conformations  with  lowest  energy,  the  difference  is 
less  than  10%.  Apparently,  the  NASA  dependence  on  confor¬ 
mation  is  relatively  small  and  is  further  narrowed  by  energy- 
delimited  restraints.  We  believe,  therefore,  that  it  is  reasonable 
to  use  the  lowest-energy  conformation  as  our  basis  for  calculat¬ 
ing  ADME  properties  in  the  current  study. 


A  QSAR  model  for  LD50  based  on  ADME  descriptors  alone 
yields  a  correlation  of  only  R2  =  0.26,  and  a  root-mean- square 
error  of  rmse  =  0.74.  This  indicates  that  substantial  discrimina¬ 
tion  among  OP  toxicants  arises  from  issues  relating  to  their 
direct  interactions  with  AChE. 

Phosphorylation.  Once  OP  compounds  access  the  active  site 
of  AChE,  they  generally  first  form  a  covalent  bond  between 
the  active  serine  by  nucleophilic  attack.  The  transient  complex 
then  forms  a  stable  tetrahedral  adduct  after  expelling  a  leaving 
group  and,  thus,  completes  the  phosphorylation  process.  Similar 
to  the  case  for  OP  metabolism,  it  is  impractical  to  do  high- 
level  quantum  chemical  transition  state  calculations  for  all  of 
the  OP  compounds  in  order  to  rigorously  predict  the  activation 
energy  barrier  and  corresponding  reaction  rate.  We  thus,  once 
again,  approximate  the  relative  propensity  of  a  given  reaction 
by  evaluating  its  enthalpy  as  a  function  of  the  initial  and  final 
states.  Note  that  the  basic  structure  of  the  AChE  active  site  prior 
to  OP  binding  does  not  vary  as  a  function  of  the  OP  itself  and 
thus  can  be  eliminated  as  a  constant  across  the  set  of  different 
OP  species.  Consequently,  for  our  thermochemical  evaluation 
of  the  reactant  state,  we  have  restricted  ourselves  to  quantum 
chemical  modeling  of  the  OP  compounds  alone,  whereas  our 
model  approximated  the  product  state  via  an  OP— Ser203 
conjugate.  In  each  case,  the  energy  corresponded  to  that  of  an 
AM  1  -optimized  geometry . 

Ligand  Aging.  The  phosphorylated  conjugate  can,  in  some 
cases,  proceed  to  the  aging  step  and  thus  complete  an  irreversible 
enzyme  inhibition.  In  general,  the  aging  involves  departure  of 
the  alkyl  group  from  the  phosphyl  alkoxy  substituent  in  the 
conjugate.  In  a  recent  study  on  nerve  agent  Tabun,  however,  it 
has  been  found  that  aging  occurs  through  the  scission  of  the 
P— N  bond  rather  than  scission  of  the  O— C  bond,  although  both 
bonds  exist  in  the  agent  (27).  Fenamiphos,  Methamidophos,  and 
Isofenphos  all  have  functional  group  similar  to  Tabun  and  thus 
should  probably  also  undergo  scission  of  the  P— N  bond  in  the 
aging  process.  The  initial  reactant  state  for  aging  is  the  Ser203— 
OP  conjugate,  and  the  product  is  the  Ser 203— OP-aged  conjugate 
plus  the  departed  aging  leaving  group,  both  of  which  have  been 
energetically  evaluated  for  each  OP  compound  again  via  AMI 
optimizations.  Considering  all  of  these  issues  in  addition  to  the 
previously  defined  metabolic  parameters  allows  us  to  derive  the 
following  relationship: 

pLD50  =  2.82  —  0.001  42£(orig)  +  0.0312£(met)  — 

(0.004  95)(NASA)  -  0.0364£(s203op)  -  0.0249£(plg)  + 
0.005  85£(s203op_aged)  -  0.0290£(alg)  (1) 

where  £(orig)  is  the  energy  of  the  parent  OP  compound,  £(met) 
is  the  energy  of  the  metabolite  compound,  NASA  is  the  negative 
accessible  surface  area  of  the  metabolite  compound,  £(s203op) 
is  the  energy  of  the  Ser203-OP  compound,  £(plg)  is  the  energy 
of  the  phosphorylation  leaving  group,  £(s203op_aged)  is  the 
energy  of  the  aged  Ser203— OP  compound,  and  E’(alg)  is  the 
energy  of  the  aging  leaving  group.  This  equation  results  in  a 
correlation  score  ( R 2  =  0.49;  rmse  =  0.62)  (Figure  3)  that  is 
improved  relative  to  models  derived  on  ADME  considerations 
alone  but  is  still  of  only  borderline  statistical  significance.  It  is 
likely  that  better  correlation  could  be  achieved  by  augmenting 
the  above  expression  with  terms  from  the  voluminous  existing 
descriptor  libraries;  however,  it  is  our  desire  to  avoid  overfitting 
the  expression  with  quantities  that  may  not  provide  obvious 
mechanistic  insight.  Indeed,  having  seven  descriptors  in  eq  1 
already  runs  some  risk  of  overfitting,  although,  mechanistically, 
it  is  difficult  to  argue  for  the  omission  of  any  of  them.  The 
intuitively  reasonable  way  of  improving  the  above  model  is  to 
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Figure  3.  Correlation  between  experimental  pLD50  and  the  prediction 
from  eq  1. 

account  for  noncovalent  precursor  (Michaelis)  binding  complex 
formation  between  the  ligand  and  enzyme  prior  to  phosphoryla¬ 
tion— a  key  to  many  inhibition  processes  and  likely  an  important 
discriminating  factor  in  the  AChE  inhibition  process. 

Subsite  Binding  Effects.  The  AChE  gorge  region  is  a 
biochemically  complex  construct  with  numerous  subsites  playing 
roles  in  the  normal  enzymatic  function  and  ultimately  also  in 
the  inhibition  process.  The  AChE  subsites  that  play  roles  in  the 
binding  of  OP  compounds  include  (a)  the  esteratic  site  contain¬ 
ing  the  active  site  serine;  (b)  the  “oxyanion  hole”  consisting  of 
residues  Glyl21,  Glyl22,  and  Ala204;  (c),  the  “anioinic  subsite” 
or  the  choline  binding  subsite,  Trp86;  (d)  the  hydrophobic  site 
for  the  alkoxy  leaving  group  of  the  substrate  containing  an 
“aromatic  patch”  that  includes  resides  Trp86,  Tyr337,  and 
Phe338;  and  (e)  the  acyl  pocket,  Phe295  and  Phe297  (28). 
Normally,  one  would  address  the  binding  affinity  by  simply 
computing  the  interaction  energy  between  OP  and  AChE. 
However,  at  the  current  stage,  (a)  the  extent  and  complexity  of 
the  AChE  active  site  render  such  calculations  impractical  at  the 
quantum  chemical  level  for  such  a  number  of  OP  compounds 
and  (b)  no  docking  programs  are  capable  of  reliably  representing 
the  formation  of  covalent  bonds.  Instead,  we  have  chosen  to 
characterize  affinity  trends  herein  via  Comparative  Molecular 
Field  Analysis  (CoMFA). 

The  most  important  step  in  CoMFA  is  to  obtain  molecular 
conformations  that  correspond  reasonably  well  to  those  observed 
for  each  ligand  in  its  bound  state.  In  previous  investigations,  it 
has  been  concluded  that  the  phosphorylation  conjugate  should 
have  the  P=0  oriented  to  the  oxyanion  hole  to  take  advantage 
of  available  hydrogen  bonding  (28).  Furthermore,  the  P— OR 
moiety  should  be  oriented  in  a  manner  that  permits  His447 
imidazolium  involvement  in  the  aging.  Such  advance  knowledge 
is  of  great  value  in  assembling  reasonable  OP  conformers  for 
the  CoMFA  and  in  interpreting  the  eventual  results.  In  the 
current  work,  all  OP  metabolites  were  aligned  to  the  crystal 
structure  of  postphosphorylation  VX.  In  the  first  step  of  the 
alignment,  the  phosphorus  center  was  aligned  with  the  VX 
phosphorus  in  1VXR  to  determine  the  anchor  point  within  the 
receptor  active  site.  Once  the  phosphorus  aligned  correctly,  the 
P  was  fixed  in  place  to  permit  rapid  alignment  of  the  rest  of 
the  molecule.  Flexible  alignment  was  done  by  weighting  the 
pharmacophore  feature  of  each  atom  with  a  Gaussian  function. 
The  aligned  structures  are  shown  in  Figure  4.  It  is  apparent  from 
the  structural  overlap  that  the  phosphorylation  leaving  group 


has  more  flexibility  than  the  other  functional  groups.  This  seems 
reasonable  given  that  the  relatively  small  space  near  the  Ser203, 
His447,  and  acyl  pocket  in  the  active  site  likely  restrains  the 
P=0  and  the  aging  leaving  group  but  does  not  impinge  on  the 
phosphorylation  leaving  group. 

To  generate  our  CoMFA  model,  an  oxygen  anion  probe  was 
used  to  calculate  the  interaction  energy  at  grid  points.  Our 
CoMFA  methodology  was  fairly  conventional  in  that  we  used 
both  electrostatic  and  van  der  Waals  interactions  to  quantify 
the  energy.  This  resulted  in  28  392  data  points  for  each  of  the 
30  compounds.  After  filtering  those  nonsignificant  terms  ac¬ 
cording  to  a  variable  importance  threshold  of  0.8,  we  were  left 
with  2586  data  points  for  each  compound  on  which  to  perform 
the  PFS  regression.  One  principle  component  was  used.  The 
resulting  correlation  between  predicted  AChE  binding  affinity 
and  experiment  pFDso  value,  shown  in  Figure  5,  gave  R 2  = 
0.892,  with  a  leave-one-out  cross-validation  score  of  g2  =  0.571. 
The  corresponding  rmse  was  0.3,  and  the  F  value  was  1 1 1.804. 
The  CoMFA  results  reveal  that  the  electrostatic  effect  is  more 
important  (fraction  of  0.657)  than  the  steric  effect  (0.343). 

To  further  analyze  the  spatial  distribution  of  favorable  and 
unfavorable  interactions  between  AChE  and  the  OP  compounds, 
the  molecular  surface  was  plotted  (Figure  6)  as  a  function  of 
all  aligned  OP  metabolites.  The  molecular  surface  shown  was 
calculated  as  the  isosurface  of  the  zero  van  der  Waals  potential. 
Regions  predicted  to  be  exposed  are  shown  in  blue,  hydrophilic 
interactions  in  red,  and  hydrophobic  in  gray.  This  molecular 
surface  is  consistent  with  the  structure  of  the  active  site  and 
the  known  OP/ AChE  inhibitive  mechanism.  The  red  region  near 
the  P=0  group  indicates  strong  negative  electrostatic  potential, 
which  corresponds  to  probable  hydrogen  bonding  within  the 
oxyanion  hole  (primarily  Glyl21  and  Glyl22).  The  field  near 
P— OR  suggests  a  possible  interaction  with  His447  in  that  the 
gray  regions  near  the  terminal  alkyl  group  indicate  favorable 
hydrophobic  interactions.  The  mixed  red,  gray,  and  blue  region 
around  the  phosphorylation  leaving  group  indicates  access  to  a 
relatively  large  space  and  dynamic  multigroup  interactions.  The 
major  force  governing  the  orientation  of  the  phosphorylation 
leaving  group  is  likely  a  hydrophobic  interaction  toward  Trp86 
and  the  entrance  of  the  gorge. 

To  further  explore  the  OP  inhibition  of  AChE,  we  also  studied 
the  interaction  between  the  phosphorylated  OP  moiety  and  the 
AChE  receptor.  The  net  electrostatic  and  van  der  Waals 
interaction  energies  were  tested  for  direct  correlation  with  acute 
toxicity;  however,  no  statistically  significant  relationship  was 
found.  As  a  result,  we  have  concluded  that  the  actual  AChE 
OP  inhibition  process  is  mainly  governed  by  the  binding  energy 
of  the  Michaelis  complex  and  the  relative  propensity  of  the 
phosphorylation  and  aging  processes. 

Unified  Model.  To  best  predict  the  acute  toxicity  of  the  OP 
compounds,  we  have  integrated  the  results  of  the  CoMFA  model 
with  those  of  the  thermochemical  and  ADME  analyses  to  arrive 
at  the  following  relationship: 

pLD50  =  -0.722  -  (1.30  x  10“4)(£(orig))  - 

1.6310“3£(met)  -  (5.00  x  10“5)(NASA)  - 
(3.60  x  10“3)(£(s203op))  +  (2.91  x  10“3)(£(plg))  + 
(5.40  x  10“4)(E(s203op_aged))  -  (4.87  x  10~3)(E(alg))  + 

(0.985)(CoMFA)  (2) 

where  the  CoMFA  term  is  the  calculated  affinity  score  from 
the  molecular  interaction  field.  Although  the  CoMFA  term 
clearly  dominates  the  overall  eq  2  and  little  improvement  is 
noted  in  the  primary  correlation  (R2  =  0.90;  rmse  =  0.301;  F 
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Figure  4.  Aligned  structures  of  the  30  OP  metabolites.  Oxygen  in  red,  carbon  in  gray,  sulfur  in  yellow,  nitrogen  in  blue,  and  hydrogen  in  white. 
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Figure  5.  Correlation  between  experimental  pLD50  and  the  predicted 
CoMFA  score. 


=  45.137;  see  Figure  7),  the  fully  integrated  model  enjoys  better 
predictivity  than  either  the  CoMFA  alone  or  the  phenomeno¬ 
logical  expression  (eq  1)  that  it  complements,  with  a  cross- 
validated  correlation  of  Q 2  =  0.82.  The  strong  cross-validated 
correlation  suggests  good  internal  consistency  that  would  not 
be  possible  with  overfitting  and,  thus,  validates  our  choice  of 
descriptors.  The  model  also  improves  substantially  on  a  previous 
study  wherein  a  correlation  ( R  =  0.7)  was  established  between 
the  acute  toxicity  (LD50)  and  in  vitro  AChE  potency  of  a  group 
of  direct-acting  OP  compounds  (29),  thus,  demonstrating  a 


unique  capacity  of  our  model  to  handle  those  OP  compounds 
that  undergo  biotransformation  prior  to  interacting  with  AChE. 

Contribution  Analysis.  To  compare  the  relative  importance 
of  the  terms  included  in  eqs  1  and  2,  we  have  relative  importance 
values  for  each  component  (i.e.,  coefficients  of  the  estimated 
normalized  linear  model)  and  have  reported  these  quantities  in 
Table  2.  The  CoMFA  term  is  clearly  the  largest  contribution  to 
eq  2,  with  the  next  largest  term,  E(plg),  having  a  relative  weight 
of  only  0.19.  Given  the  fact  that  the  CoMFA  predictions,  taken 
by  themselves,  correlate  strongly  ( R  =  0.71)  with  eq  1  (the 
phenomenological  QSAR  developed  without  a  CoMFA  contri¬ 
bution)  it  appears  that  much  of  this  perceived  dominance  must 
actually  result  from  functional  overlap  between  features  repre¬ 
sented  in  the  CoMFA  model  and  those  inherent  in  other 
descriptors.  Adding  a  CoMFA  term  to  the  QSAR  thus  introduces 
some  redundancy  among  descriptors  and  thereby  reduces  the 
proportional  weights  of  any  other  descriptor  that  exhibits  some 
degree  of  linear  dependence  with  the  CoMFA  (in  this  case  all 
of  them).  The  large  reduction  in  weight  of  the  NASA  term,  for 
example,  implies  that  it  is  so  nearly  completely  described  within 
the  CoMFA  as  to  contribute  negligibly  to  the  expression.  This 
is  not  especially  surprising,  since  CoMFA  models  frequently 
provide  very  good  representations  of  three-dimensional  elec¬ 
trostatic  distributions.  While  we  retain  the  NASA  term  within 
eq  2  for  the  sake  of  direct  comparison  to  eq  1,  we  note  that 
retraining  the  other  descriptors  in  the  absence  of  NASA  yields 
an  expression  with  values  for  R 2,  Q 2,  and  coefficients  for  the 
other  descriptors  that  are  identical  (within  the  significant  fig¬ 
ures  reported  herein)  to  those  reported  for  eq  2.  Furthermore, 
since  descriptors  based  on  molecular  energy  simultaneously 
represent  multiple  terms  that  may  have  varying  levels  of 
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Figure  6.  Interaction  of  molecular  surface  for  OP  compounds.  Red  regions  indicate  hydrophilic  interactions,  blue  regions  indicate  exposed  regions, 
and  gray  regions  experience  hydrophobic  interactions. 
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Figure  7.  Correlation  between  the  experimental  and  predicted  pLD50 
for  the  30  OP  compounds  for  the  fully  integrated  model  (eq  2). 

Table  2.  Relative  Importance  of  the  Descriptor  in  Eqs  1  and  2 

terms 

relative  weights 

eq  1 

eq  2 

change 

F(plg) 

0.72 

0.19 

-74% 

F(alg) 

0.22 

0.08 

-64% 

£(s203op) 

0.38 

0.08 

-79% 

F(met) 

1.0 

0.12 

-88% 

NASA 

0.17 

0.004 

-98% 

F(s203op_aged) 

0.05 

0.01 

-80% 

F(orig) 

0.05 

0.01 

-80% 

CoMFA 

- 

1.00 

+  100% 

correlative  redundancy  relative  to  the  CoMFA  (i.e.,  molecular 
electrostatic  and  van  der  Waals  effects  are  also  represented  in 
CoMFA,  but  bond  energy  is  not),  the  coefficient  for  that 


descriptor  (its  coefficient  weight  and  possibly  even  sign)  within 
a  QSAR  may  change  upon  integration  of  CoMFA  into  the 
model.  Because  of  such  effects,  eq  2,  although  a  predictive 
tool  of  exceptional  value,  is  less  amenable  to  intuitive  physi¬ 
cal  interpretation  than  the  purely  phenomenological  eq  1.  We 
will  thus  restrict  our  interpretive  analysis  to  trends  observed  in 
eq  1. 

In  eq  1,  £(met)  and  £(plg)  are  assigned  the  largest  weights, 
while  £(s203op_aged)  and  £(orig)  appear  to  be  the  least  impor¬ 
tant  descriptors.  This  observation  is  consistent  with  the  mech¬ 
anism  of  actual  OP  toxicity.  For  £(plg)  versus  £(s203op_aged), 
the  physical  rationale  for  these  weights  is  quite  clear:  structures 
of  the  aged  complexes  vary  fairly  little  across  the  set  of  30 
toxins,  whereas  the  phosphorylation  leaving  group  contains 
much  of  the  molecular  diversity  present  in  this  collection.  For 
£(met)  versus  £(orig),  an  important  kinetic  relationship  is 
suggested:  if  the  parent  compound  energies  ultimately  have 
relatively  small  contributions  on  the  QSAR,  it  suggests  that  their 
respective  metabolic  reactions  may  all  share  a  similar  activation 
energy,  whereas  the  importance  of  i^met)  suggests  that  the 
reverse  activation  barrier  may  have  substantially  greater  diver¬ 
sity.  £(met)  also  must  be  considered  as  descriptor  of  relevance 
to  the  subsequent  phosphorylation  process  (i.e.,  in  this  case 
modeling  the  reactant);  thus  in  effect,  it  contributes  twice  to 
the  equation. 


Conclusions 

Molecular  modeling  has  been  used  herein  to  investigate  the 
mechanism  of  acute  toxicity  of  the  OP  compounds  via  the 
construction  of  a  number  of  QSAR  models  to  account  for 
available  LD50  data.  A  purely  phenomenological  QSAR  model 
describing  the  toxic  process  as  a  function  of  ADME  and  covalent 
bonding  processes  was  constructed  as  a  seven-parameter  model 
wherein  ADME  effects  were  modeled  as  the  reactant  and 
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product  energies  for  the  incumbent  biotransformation  reaction, 
plus  the  ligand  negative  accessible  surface  area,  and  the  covalent 
phosphorylation  and  aging  reactions  were  approximated  by 
thermochemical  characterization.  This  model  achieved  a  cor¬ 
relation  of  R 2  =  0.49  relative  to  experiment,  suggesting  that 
the  relationship  is  statistically  valid  but,  lacking  an  account  of 
the  initial  ligand— receptor  Michaelis  complex  formation,  likely 
not  of  great  predictive  value. 

Given  the  unreliability  of  using  docking  methods  to  evaluate 
affinities  for  precovalent  Michaelis  complexes,  we  chose  to 
account  for  this  effect  via  a  CoMFA  pseudointeraction  field. 
Without  suitable  in  vitro  affinity  data  covering  our  full  training 
set,  we  chose  to  fit  the  CoMFA  model  directly  to  the  LD50, 
achieving  a  surprisingly  strong  correlation  with  experiment  (. R 2 
=  0.89),  thus,  indicating  that  the  pharmacophoric  fit  between 
receptor  and  ligand  is  a  key  factor  in  determining  toxic  inhibition 
of  the  AChE  enzyme.  However  the  CoMFA' s  marginal  predic- 
tivity  ( Q 2  =  0.57)  clearly  reflects  a  need  for  a  more  rigorous 
account  of  ADME  effects  and  subsequent  covalent  binding 
processes.  With  all  parameters  considered  together  (ADME, 
covalent,  and  CoMFA),  both  correlation  (. R 2  =  0.90)  and 
predictivity  ( Q 2  =  0.82)  are  strong. 

While  the  unified  model  is  of  substantial  predictive  value,  it 
does  come  at  the  expense  of  interpretive  ease:  while  spatial 
aspects  of  the  CoMFA  reproduce  AChE  inhibition  phenomena 
that  have  been  previously  observed  from  structure-based 
analysis,  correlative  redundancy  between  the  CoMFA  and  the 
various  phenomenological  descriptors  tends  to  produce  nonin- 
tuitive  weights  and  signs  that  partially  obscure  the  role  of  the 
various  terms.  For  this  reason,  the  purely  phenomenological 
(non-CoMFA)  model  remains  of  some  rhetorical  value,  provid¬ 
ing  some  guidance  as  to  the  relative  importance  of  different 
ADME  and  covalent  binding  effects. 

The  fact  that  both  our  phenomenological  and  composite 
models  exhibit  reasonable  correlation  with  existing  LD50  data 
despite  approximating  the  efficacy  of  multiple  enzymatic  and 
covalent  inhibition  reactions  via  thermochemistry  rather  than 
kinetic  activation  barriers  may  suggest  that  the  reactive  processes 
must  largely  achieve  equilibrium  within  the  time  frame  under 
which  the  in  vivo  studies  were  carried  out— a  fact  that  puts  more 
weight  on  the  relative  energetics  of  product  and  reaction  than 
on  the  actual  reaction  barrier.  This  assessment  deserves  further 
analysis  through  explicit  quantum  chemical  transition  state 
calculations  that  are  planned  as  a  future  study.  Our  fairly 
simplistic  ADME  model  is  another  area  for  potential  improve¬ 
ment.  Indications  are  that  our  metabolic  representation  is 
reasonable;  however,  the  blanket  representation  of  other  absorp¬ 
tion,  distribution,  and  excretion  effects  via  single  NASA  is  a 
major  approximation.  While  the  CoMFA  appeared  to  be 
surprisingly  adept  at  encoding  such  ADME  trends,  it  is  highly 
unlikely  that  it  recovered  all  relevant  effects.  For  example,  one 
factor  of  significant  possible  interest  to  AChE  inhibition  research 
that  may  not  have  been  adequately  accounted  for  is  blood— 
brain  barrier  (BBB)  transmission.  Identification  of  a  convenient 
and  apt  BBB  descriptor  may  further  improve  the  model. 

Rigorously  addressing  all  potential  model  refinements  may 
constitute  a  computationally  unwieldy  paradigm.  While  such 
refinement  should  nonetheless  be  pursued  for  the  sake  of  a 
complete  understanding  of  the  OP  toxicity  mechanism,  it  is  our 
short-term  hope  that  the  relative  computational  simplicity  of 
this  current  model,  as  well  as  its  reasonable  predictive  capacity, 
will  permit  rapid  and  reliable  evaluation  of  other  prospective 
toxins.  The  simpler  phenomenological  model  also  pinpoints 
trends  in  the  relative  importance  of  different  factors  on  the  acute 


toxicity  but  does  make  substantial  approximations  and  leaves 
some  details  of  the  specific  effects  rather  opaque.  In  the  interests 
of  full  understanding  of  the  acute  toxicity,  we  are  considering 
some  additional  factors,  including  the  influence  of  water  on  the 
catalytic  reaction  within  AChE,  dynamic  conformational  changes 
in  the  active  site,  and  transport  kinetics  governing  how  the  OP 
compounds  access  the  catalytic  triad.  In  addition,  mutagenesis 
work  suggests  that  the  peripheral  anionic  site  may  affect  the 
rate  constant  by  interacting  with  the  cationic  aging  leaving  group 
at  Asp74  in  human  AChE  (30).  Finally,  the  details  relating  to 
the  binding  of  the  anionic  phosphorylation  leaving  group  may 
also  prove  to  be  of  interest. 

Acknowledgment.  This  work  was  supported  by  an  Army 
Research  Office  Short  Term  Innovative  Research  award  (No. 
W911NF-05-01-0181).  The  authors  thank  the  Department  of 
Defense  High  Performance  Computing  Modernization  Program 
(Project  No.ARLAP00583C91)  and  the  National  Center  for 
Supercomputing  Applications  (Project  No.  MCB030011N)  for 
computational  resources. 

Supporting  Information  Available:  All  terms  used  in  eqs  1 
and  2  for  each  compound  have  been  provided  as  Supporting 
Information.  All  structures  for  the  toxins,  metabolites,  complexes, 
and  leaving  groups  are  available  upon  request.  This  material  is 
available  free  of  charge  via  the  Internet  at  http://pubs.acs.org. 

References 

(1)  Flanagan,  J.,  and  Jones,  A.  L.  (2001)  Antidotes,  pp  87—114,  Taylor 
and  Francis,  London  and  New  York. 

(2)  Ehrich,  M.  (1998)  Organophosphates.  In  Encyclopedia  of  Toxicology 
(Wexler,  P.,  Ed.)  pp  467—471,  Academic  Press,  San  Diego,  CA. 

(3)  Sidell,  F.  R.,  and  Borak,  J.  (1992)  Chemical  warfare  agents:  II.  Nerve 
agents.  Ann.  Emerg.  Med.  21,  865—871. 

(4)  Marrs,  T.  C.  (1993)  Organophosphate  poisoning.  Pharmacol.  Ther. 
58,  51-66. 

(5)  Sussman,  J.  L.,  Harel,  M.,  Frolow,  F.,  Oefner,  C.,  Goldman,  A.,  Toker, 
L.,  and  Silman,  I.  (1991)  Atomic  structure  of  acetylcholinesterase  from 
Torpedo  califomica:  a  prototypic  acetylcholine-binding  protein. 
Science  253,  872—879. 

(6)  Kovarik,  Z.,  Radic,  Z.,  Berman,  H.  A.,  Simeon-Rudolf,  V.,  Reiner, 
E.,  and  Taylor,  P.  (2003)  Acetylcholinesterase  active  centre  and  gorge 
conformations  analysed  by  combinatorial  mutations  and  enantiomeric 
phosphonates.  Biochem.  J.  373,  33—40. 

(7)  Gilson,  M.  K.,  Straatsma,  T.  P.,  McCammon,  J.  A.,  Ripoll,  D.  R., 
Faerman,  C.  H.,  Axelsen,  P.  H.,  Silman,  I.,  and  Sussman,  J.  L.  (1994) 
Open  “back  door”  in  a  molecular  dynamics  simulation  of  acetylcho¬ 
linesterase.  Science  263,  1276—1278. 

(8)  Pope,  C.  N.  (1999)  Organophosphorus  pesticides:  do  they  all  have 
the  same  mechanism  of  toxicity?  1.  Toxicol.  Environ.  Health,  Part  B 
2,  161-181. 

(9)  Knaak,  J.  B.,  Dary,  C.  C.,  Power,  F.,  Thompson,  C.  B.,  and  Blancato, 
J.  N.  (2004)  Physicochemical  and  biological  data  for  the  development 
of  predictive  organophosphorus  pesticide  QSARs  and  PBPK/PD 
models  for  human  risk  assessment.  Crit.  Rev.  Toxicol.  34,  143—207. 

(10)  Barril,  X.,  Orozco,  M.,  and  Luque,  F.  J.  (2001)  Towards  improved 
acetylcholinesterase  inhibitors:  a  structural  and  computational  ap¬ 
proach.  Mini  Rev.  Med.  Chem.  1,  255—266. 

(11)  Kua,  J.,  Zhang,  Y.,  and  McCammon,  J.  A.  (2002)  Studying  enzyme 
binding  specificity  in  acetylcholinesterase  using  a  combined  molecular 
dynamics  and  multiple  docking  approach.  J.  Am.  Chem.  Soc.  124, 
8260-8267. 

(12)  Guo,  J.-X.,  Hurley,  M.  M.,  Wright,  J.  B.,  and  Lushington,  G.  H.  (2004) 
A  docking  score  function  for  estimating  ligand-protein  interactions: 
application  to  acetylcholinesterase  inhibition.  J.  Med.  Chem.  47,  5492— 
5500. 

(13)  Berman,  H.  M.,  Westbrook,  J.,  Feng,  Z.,  Gilliland,  G.,  Bhat,  T.  N., 
Weissig,  H.,  Shindyalov,  I.  N.,  and  Bourne,  P.  E.  (2000)  The  protein 
data  bank.  Nucleic  Acids  Res.  28,  235—242. 

(14)  Millard,  C.  B.,  Koellner,  G.,  Ordentlich,  A.,  Shafferman,  A.,  Silman, 
I.,  and  Sussman,  J.  L.  (1999)  Reaction  products  of  acetylcholinesterase 
and  VX  reveal  a  mobile  histidine  in  the  catalytic  triad.  J.  Am.  Chem. 
Soc.  121,  9883-9884. 


216  Chem.  Res.  Toxicol,  Vol.  19,  No.  2,  2006 


Guo  et  al. 


(15)  Bourne,  Y.,  Taylor,  P.,  Radic,  Z.,  and  Marchot,  P.  (2003)  Structural 
insights  into  ligand  interactions  at  the  acetylcholinesterase  peripheral 
anionic  site.  EMBO  J.  22,  1  —  12. 

(16)  Legay,  C.,  Bon,  S.,  Vernier,  P.,  Coussen,  F.,  and  Massoulie,  J.  (1993) 
Cloning  and  expression  of  a  rat  acetylcholinesterase  subunit:  genera¬ 
tion  of  multiple  molecular  forms  and  complementarity  with  a  Torpedo 
collagenic  subunit.  J.  Neurochem.  60,  337—346. 

(17)  Weiner,  S.  J.,  Kollman,  P.  A.,  Nguyen,  D.  T.,  and  Case,  D.  A.  (1986) 
An  all  atom  force  field  for  simulations  of  proteins  and  nucleic  acids. 
J.  Comput.  Chem.  7,  230—252. 

(18)  Halgren,  T.  A.  (1999)  MMFF  VI.  MMFF94s  option  for  energy 
minimization  studies.  J.  Comput.  Chem.  20,  720—729. 

(19)  Molecular  Opperating  Enviroment  (MOE),  version  2004.03,  Chemical 
Computing  Group,  Inc.,  Montreal,  Quebec,  Canada,  Mar,  2004. 

(20)  Dewar,  M.  J.  S.,  Zoebisch,  E.  G.,  Healy,  E.  F.,  and  Stewart,  J.  J.  P. 
(1995)  AMI:  a  new  general  purpose  quantum  mechanical  model.  J. 
Am.  Chem.  Soc.  107,  3902—3909. 

(21)  Stewart,  J.  J.  P.  (1993)  MOPAC7,  version  7,  Stewart  Computational 
Chemistry,  Colorado  Springs,  CO. 

(22)  Ferguson,  D.  M.,  and  Raber,  D.  J.  (1989)  A  new  approach  to  probing 
conformational  space  with  molecular  mechanics:  random  incremental 
pulase  search.  J.  Am.  Chem.  Soc.  Ill,  4371—4378. 

(23)  SIMCA-P,  Umereics  AB,  Umea,  Sweden,  July,  2001. 

(24)  Van  der  Waterbeemed,  H.,  and  Kansy,  M.  (1992)  Hydrogen-bonding 
capacity  and  brain  pentration.  Chimia  46,  299—303. 

(25)  Hou,  T.  J.,  Zhang,  W.,  Xia,  K.,  Qiao,  X.  B.,  and  Xu,  X.  J.  (2004) 
ADME  evaluation  in  drug  discovery.  5.  Correlation  of  Caco-2 
permeation  with  simple  molecular  properties.  J.  Chem.  Inf.  Comput. 
Sci.  44,  1585-1600. 

(26)  Liu,  R.,  Sun,  H.,  and  So,  S.-S.  (2001)  Development  of  quantitative 
structure— property  relationship  models  for  early  ADME  evaluation 
in  drug  discovery.  2.  Blood— brain  barrier  penetration.  J.  Chem.  Inf. 
Comput.  Sci.  41,  1623—1632. 

(27)  Barak,  D.,  Ordentlich,  A.,  Kaplan,  D.,  Barak,  R.,  Mizrahi,  D., 
Kronman,  C.,  Segall,  Y.,  Velan,  B.,  and  Shafferman,  A.  (2000) 
Evidence  for  P— N  bond  scission  in  phosphoroamidate  nerve  agent 
adducts  of  human  acetylcholinesterase.  Biochemistry  39,  1156—1161. 

(28)  Ordentlich,  A.,  Barak,  D.,  Kronman,  C.,  Benschop,  H.  P.,  De  Jong, 
L.  P.,  Ariel,  N.,  Barak,  R.,  Segall,  Y.,  Velan,  B.,  and  Shafferman,  A. 
(1999)  Exploring  the  active  center  of  human  acetylcholinesterase  with 
stereomers  of  an  organophosphorus  inhibitor  with  two  chiral  centers. 
Biochemistry  38,  3055—3066. 

(29)  Holmstedt,  B.  (1963)  Structure— activity  relationship  of  the  organo¬ 
phosphorus  anticholinesterase  agents.  Handb.  Exp.  Toxicol.  15,  428— 
485. 

(30)  Ordentlich,  A.,  Barak,  D.,  Sod-Moriah,  G.,  Kaplan,  D.,  Mizrahi,  D., 
Segall,  Y.,  Kronman,  C.,  Karton,  Y.,  Lazar,  A.,  Marcus,  D.,  Velan, 
B.,  and  Shafferman,  A.  (2004)  Stereoselectivity  toward  VX  is 
determined  by  interactions  with  residues  of  the  acyl  pocket  as  well  as 
of  the  peripheral  anionic  site  of  AChE.  Biochemistry  43,  1 1255 — 
11265. 

(31)  Yamanaka,  S.,  Ohta,  K.,  Tamita,  Y.,  Takayanagi,  A.,  Nomura,  T.,  and 
Takaesu,  Y.  (1996)  Effects  on  acute  organophosphorus  poisoning  in 
rats  in  aging  and  solubility  of  organophosphates.  Environ.  Health  Prev. 
Med.  1,  119-127. 

(32)  Lin,  S.  N.,  Chen,  C.  Y.,  Murphy,  S.  D.,  and  Caprioli,  R.  M.  (1980) 
Quantitative  high-performance  liquid  chromatography  and  mass 
spectrometry  for  the  analysis  of  the  in  vitro  metabolism  of  the 
insecticide  azinphos-methyl  (guthion)  by  rat  liver  homogenates.  J. 
Agric.  Food  Chem.  28,  85—88. 

(33)  Lucier,  G.  W.,  and  Menzer,  R.  E.  (1970)  Nature  of  oxidative 
metabolites  of  dimethoate  formed  in  rats,  liver  microsomes,  and  bean 
plants.  J.  Agric.  Food  Chem.  18,  698—704. 

(34)  Ivey,  M.  C.,  and  Mann,  H.  D.  (1975)  Gas— liquid  chromatographic 
determination  of  ethion,  ethion  monooxon,  and  ethion  dioxon  in  tissues 
of  turkeys  and  cattle.  J.  Agric.  Food  Chem.  23,  319—321. 


(35)  Krueger,  H.  R.,  and  O’Brien,  R.  D.  (1959)  Relation  between 
metabolism  and  differential  toxicity  of  malathion  in  insects  and  mice. 
J.  Econ.  Entomol.  52,  1063—1067. 

(36)  Chopade,  H.  M.,  and  Dauterman,  W.  C.  (1981)  Studies  on  the  in  vitro 
metabolism  of  methidathion  by  rat  and  mouse  liver.  Pestic.  Biochem. 
Physiol.  15,  105-119. 

(37)  Venkatesh,  K.,  Levi,  P.  E.,  Inman,  A.  O.,  Monteiro-Riveiere,  N.  A., 
Misra,  R.,  and  Hodgson,  E.  (1992)  Enzymic  and  immunohistochemical 
studies  on  the  role  of  cytochrome  p450  and  the  flavin-containing 
monooxygenase  of  mouse  skin  in  the  metabolism  of  pesticides  and 
other  xenobiotics.  Pestic.  Biochem.  Physiol.  43,  53—66. 

(38)  De  Potas,  G.  M.,  and  de  D’Angelo,  A.  M.  P.  (1993)  Phosphoinositide 
phosphorylation  and  shape  changes  produced  by  phosmet-oxon  in 
human  erythrocytes.  Comp.  Biochem.  Physiol.,  Part  C:  Pharmacol., 
Toxicol.  Endocrinol.  106,  561—566. 

(39)  Clark,  D.  E.,  Ivie,  G.  W.,  Crookshank,  H.  R.,  and  Devaney,  J.  A. 
(1979)  Effects  of  sulprofos  and  its  sulfoxide  and  sulfone  metabolites 
on  laying  hens  fed  the  compounds  in  the  diet.  J.  Agric.  Food  Chem. 
27,  103-107. 

(40)  Li,  J.-T.,  Sheng,  S.-J.,  and  Du,  X.-L.  (1999)  Metabolism  of  terbufos 
in  rat  liver.  J.  Occup.  Health  41,  62—68. 

(41)  Fukuto,  T.  R.  (1978)  Insecticide  metabolism  and  mode  of  action.  Pure 
Appl.  Chem.  50,  9—10. 

(42)  Beike,  J.,  Ortmann,  C.,  Meiners,  T.,  Brinkmann,  B.,  and  Kohler,  H. 
(2002)  LC_MS  determination  of  oxydemeton-methyl  and  its  main 
metabolite  demeton-s-methylsulfon  in  biological  sepcimens-applicatin 
to  a  forensic  case.  J.  Anal.  Toxicol.  25,  308—312. 

(43)  Gotoh,  M.,  Sakata,  M.,  Endo,  T.,  Hayashi,  H.,  Seno,  H.,  and  Suzuki, 
O.  (2001)  Profenofos  metabolites  in  human  poisoning.  Forensic  Sci. 
Int.  116,  221-226. 

(44)  Abu-Qare,  A.  W.,  and  Abou-Donia,  M.  B.  (2001)  Quantification  of 
nicotine,  chlorpyrifos  and  their  metabolites  in  rat  plasma  and  urine 
using  high-performance  liquid  chromatography.  J.  Chromatogr.,  B: 
Biomed.  Sci.  Appl.  757,  295  —  300. 

(45)  Osman,  A.  Z.,  Zayed,  S.  M.  A.  F.  D.,  and  Hazzaa,  N.  I.  (1986)  Fate 
and  metabolism  of  the  radiolabeled  insecticide  coumaphos  in  egyptian 
lactating  goats.  Isot.  Radiat.  Res.  18,  139—146. 

(46)  Kappers,  W.  A.,  Edwards,  R.  J.,  Murray,  S.,  and  Boobis,  A.  R.  (2001) 
Diazinon  is  activated  by  CYP2C19  in  human  liver.  Toxicol.  Appl. 
Pharmacol.  177,  68—76. 

(47)  Sultatos,  L.  G.  (1991)  Metabolic  activation  of  the  organophosphorus 
insecticides  chlorpyrifos  and  fenitrothion  by  perfused  rat  liver. 
Toxicology  68,  1—9. 

(48)  Cabras,  P.,  Plumitallo,  A.,  and  Spanedda,  L.  (1991)  High-performance 
liquid  chromatographic  separation  of  fenthion  and  its  metabolites.  J. 
Chromatogr.  540,  406—410. 

(49)  Bakke,  J.  E.,  and  Price,  C.  E.  (1976)  Metabolism  of  O,  O-dimethyl- 
O— (3,5,6-trichloro-2-pyridyl)phosphorothioate  in  sheep  and  rats  and 
of  3,5,6-trichloro-2-pyridinol  in  sheep.  J.  Environ.  Sci.  Health,  Part 
B  11,  9-22. 

(50)  Abu-Qare,  A.  W.,  and  Abou-Donia,  M.  B.  (2000)  Urinary  excretion 
of  metabolites  following  a  single  dermal  dose  of  [14C]  methyl  parathion 
in  pregnant  rats.  Toxicology  150,  119—127. 

(51)  Soranno,  T.  M.,  and  Sultatos,  L.  G.  (1992)  Biotransformation  of  the 
insecticide  parathion  by  mouse  brain.  Toxicol.  Lett.  60,  27—37. 

(52)  Berealey,  C.  J.,  and  Lawrence,  D.  K.  (1979)  High-performance  liquid 
chromatography  of  pirimiphos  methyl  and  five  metabolites.  J.  Chro¬ 
matogr.  168,  461—469. 

(53)  Gorder,  G.  W.,  Kirino,  O.,  Hirashima,  A.,  and  Casida,  J.  E.  (1986) 
Bioactivation  of  isofenphos  and  analogs  by  oxidative  N-dealkylation 
and  desulfuration.  J.  Agric.  Food  Chem.  34,  941—947. 

TX050090R 


